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The  numerical  and  analytical  study  of  bifurcation  and  multicellular 
flow  instability  due  to  natural  convection  between  narrow  horizontal 
isothermal  cylindrical  annuli  at  high  Rayleigh  numbers 


Daniel  Bartholemew  Fant 

Under  the  supervision  of  Joseph  M.  Prusa 
From  the  Department  of  Mechanical  Engineering 
Iowa  State  University 

This  resea rch^ef fort  deals  with  a  numerical  and  analytical  study 
of  multicellular  flow  instability  due  to  natural  convection  between 
narrow  horizontal  isothermal  cylindrical  annuli. 

Buoyancy- induced  steady  or  unsteady  flow  fields  between  the  annuli 
are  determined  using  the  Boussinesq  approximated  two-dimensional  {2-0) 
Navier-Stokes  equations  and  the  viscous-dissipation  neglected  thermal- 
energy  equation.  The  vorticity^stream  function  formulation  of  the 
Navier-Stokes  equations  is  adopted. 

Both  thermal  and  hydrodynamic  instabilities  are  explored.  An 
asymptotic  expansion  theory  is  applied  to  the  Navier-Stokes  equations  in 
the  double-limit  of  Rayleigh  number  approaching  infinity  and  gap  width 
approaching  zero.  This  double-limiting  condition  reduces  the  governing 
equations  to  a  set  of  Cartesian-1  ike  boundary-layer  equations.  These 
equations  are  further  simplified  by  considering  the  extreme  limits  of 
Pr  •*  «  and  Pr  -*■  0.  The  former  limit  yields  an  energy  equation  which 
retains  the  nonlinear  convective  terms,  while  the  vorticity  equation 
reduces  to  a  Stokes-flow  equation,  signifying  the  potential  for  thermal 
instability.  In  the  latter  limit,  the  nonlinear  terms  in  the  vorticity 


equation  remain,  while  the  energy  equation  collapses  to  a  one-dimensional 
conduction  equation,  signifying  the  potential  for  hydrodynamic 
instability. 

Thermal  instability  of  air  near  the  top  portions  of  narrow  annuli 
is  considered  for  various  size  small  gap  widths.  For  these  narrow  gaps, 
the  Rayleigh  numbers  corresponding  to  the  onset  of  steady  multicellular 
flow  are  predicted.  Numerical  solutions  of  the  2-D  Navier-Stokes 
equations  also  yield  hysteresis  behavior  for  the  two-to-six  and  two-to- 
four  cellular  states,  with  respect  to  diameter  ratios  of  1.100  and  1.200. 
In  contrast,  an  unsteady  hydrodynamic  multicellular  instability  is 
experienced  near  the  vertical  sections  of  narrow  annuli  when  the  Pr  -*■  0 
boundary-layer  equations  are  solved  numerically. 

In  addition,  analytical  steady-state  perturbative  solutions  to  the 
boundary- layer  equations  are  obtained.  These  results  compare  favorably 
to  related  numerical  solutions  of  both  the  Navier-Stokes  and  the  Pr  -*■  0 
simplified  equations. 

In  all  cases,  finite-differenced  solutions  to  the  governing 
equations  are  obtained  using  a  stable  second-order,  fully-implicit 
time-accurate  Gauss-Seidel  iterative  procedure. 
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1.  INTRODUCTION 


From  a  practical  standpoint,  the  study  of  natural  convec*-  - 
between  horizontal  isothermal  concentric  cylinders  has  a  wise  .a-  . 

of  technological  applications,  ranging  from  nuclear  reactors  arc 
thermal  storage  systems  to  cooling  of  electronic  components,  a1*- 
fuselage  insulation,  underground  electrical  transmission  lines  arc  eve' 
the  flow  in  the  cooling  passages  of  turbine  blades  (Tsui  and  Tremblay, 
1984). 

However,  in  a  different  perspective,  the  work  set  forth  in  this 
research  effort  was  undertaken  to  gain  a  more  practical  understanding 
of  the  effects  of  nonlinearity  with  regard  to  natural  convective  flow 
instabilities.  Still  not  well  understood  is  the  influence  of  °randtl 
number  variations  on  the  nonlinear  processes  involved  in  triggering 
either  thermal  or  hydrodynamic  types  of  instabilities.  Also  of 
relevance  is  the  aspect  of  nonuniqueness,  which  allows  the  possibility 
of  hysteresis  behavior  associated  with  thermal -convective 
instabilities.  These  important  issues  and  concerns  are  addressed  and 
studied  in  this  thesis. 

In  light  of  the  above,  the  work  in  this  thesis  is  separated  into 
three  main  areas. 

First,  a  stable  second-order  finite-difference  solution  to  the 
2-D  Navier-Stokes  equations  is  implemented  in  order  to  investigate 
possible  hysteresis  behavior  in  relation  to  multicellular  thermal 
instability  near  the  top  of  narrow  horizontal  annuli  for  air. 


The  second  area  of  work  involves  an  asymptotic  expansion  theory 
applied  to  the  2-D  Navier-Stokes  equations  in  the  double-limit  of 
Rayleigh  number  approaching  infinity  and  gap  width  approaching  zero. 
In  this  double-limit,  the  Navier-Stokes  equations  are  reduced  to 
Cartesian-like  boundary-layer  equations.  Analytical  steady-state 
solutions  to  these  simplified  equations  are  also  obtained,  and  the 
results  are  compared  to  related  2-D  Navier-Stokes  numerical  data. 
Moreover,  in  order  to  obtain  further  insight  into  the  nonlinearity 
associated  with  extreme  Prandtl  number  variations,  limiting  boundary- 
layer  equations  for  Pr  ->  0  and  Pr  «  are  derived. 

Thirdly,  the  Pr  ^  0  simplified  boundary-layer  equations  are 
solved  numerically  to  investigate  the  full  effects  of  nonl inearity, 
which  were  believed  to  cause  an  unsteady  hydrodynamic  multicellular 
instability  between  the  vertical  portions  of  narrow  horizontal  annuli 

Before  delving  into  deep  analysis  and  discussion  of  these  topics 
an  intensive  literature  review  and  derivation  of  the  applicable 
governing  equations  are  presented  in  Chapters  2  and  3,  respectively. 
Finally,  the  key  results  and  conclusions  of  this  work  are  given  in 
Chapters  6  and  7. 
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2.  LITERATURE  REVIEW 

Natural  convection  phenomena  between  horizontal  isothermal 
concentric  cylinders  have  been  scrutinized  experimentally,  analytically, 
and  numerically  throughout  the  past  few  decades.  In  the  1960s,  most 
work  was  experimental  in  nature.  In  the  '70s  and  ‘80s,  many  numerical 
studies  dominated  the  literature  due  to  the  advent  of  the  modern  computer. 
Analytical  studies,  in  general,  have  been  much  more  limited. 

This  thesis  concentrates  on  the  high  Rayleigh  number/smal 1 -gap 
flow  regime.  It  has  been  found  that  while  numerical  studies  in  the  low 
to  moderate  Rayleigh  number  range  are  quite  abundant  and  agree  rather 
favorably  with  related  experimental  work,  the  numerical  work  in  the 
high  Rayleigh  number  multicellular  flow  regime  (usually  associated  with 
narrow  gaps)  has  been  much  less  exhaustive.  This  is  due  to  the  fact 
that  many  computational  schemes  either  could  not  resolve  the  transition 
to  the  laminar  multicellular  flow  field,  or  they  became  unstable  just 
prior  to  it.  Analytical  approaches,  especially  with  regard  to  the 
high  Rayleigh  number/smal 1 -gap  flow  regime,  have  been  virtually 
unexplored. 

To  better  understand  how  the  various  studies  were  conducted  and 
evolved,  this  literature  survey  will  be  essentially  divided  into  three 
main  categories;  namely,  the  analytical  analyses,  the  experimental 
approaches,  and  the  numerical  studies.  In  addition,  three  sections  will 
be  included  near  the  end  of  this  chapter  to  further  support  some  of  the 
assumptions  and  findings  relevant  to  this  present  study.  These 
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sections  will  touch  upon  variable  fluid  property  effects,  flow 
bifurcation,  and  natural  convective  flow  between  vertical  slots. 

Although  many  of  the  reviews  will  relate  to  the  pretransition,  non- 
multicellular  flow  studies,  their  inclusion  is  necessary  to  fully 
appreciate  some  of  the  follow-up  work  described  in  this  thesis. 

The  basic  flow  field  normally  encountered  between  horizontal 
isothermal  concentric  cylinders  (see  Figure  6.3a)  is  the  bicellular 
kidney-shaped  pattern  that  results  strictly  from  buoyancy  effects: 
density  differences  occur  due  to  the  inner  cylinder  being  hotter  than 
the  outer  one.  The  lighter,  hotter  fluid  begins  to  rise  in  a  boundary- 
layer  manner  near  the  warmer  inner  cylinder  while  ascending  more 
uniformly  in  the  inviscid  core  region  near  the  center  of  the  annulus. 
Finally,  it  separates  and  impinges  upon  the  top  of  the  outer  cylinder 
via  the  thermal  plume.  The  cooler,  more  dense  fluid,  then,  descends 
along  the  outer  cylinder  and  regains  its  upward  cyclic  ascent  near  the 
lower  portion  of  the  annulus.  A  similar  kidney-shaped  flow  pattern  will 
also  occur  when  the  inner  cylinder  is  cooled  and  the  outer  one  heated. 
Other  interesting  flow  patterns  are  possible  and  will  be  discussed  in 
the  sections  that  follow.  Symbols  will  be  either  directly  defined  or 
their  meaning  may  be  found  in  the  Nomenclature. 

2.1.  Analytical  Studies 

Eight  analytical  studies  pertaining  to  natural  convective  flow 
between  horizontal  concentric  cylinders  were  found  in  the  literature 
search.  Of  the  eight,  seven  involved  perturbation  methods  and  only 


one  dealt  with  linearized  stability  theory.  All  but  one  related  to  the 
relatively  small  Rayleigh  number  domain. 

Mack  and  Bishop  (1968)  used  a  Rayleigh  number  power  series 
expansion  to  obtain  a  steady-state  solution  for  natural  convection 
between  2-D  horizontal  isothermal  cylinders.  The  dimensionless 
vorticity-transport  and  energy  equations  were  assumed  to  govern  the  flow 
field  in  this  analysis.  They  expanded  temperature  (T)  and  the  stream 
function  (40  in  the  following  nweier: 

oo 

T  =  Z  AJT.(r,6)  (2.1a) 

j=0  3 


t  =  Z  A3V(r,0)  (2.1b) 

j=l  3 

where  A  signified  the  Rayleigh  number  and  0  the  angular  coordinate. 

The  first  term  in  their  expansion  represented  the  creeping-flow 
solution  and  was  in  agreement  with  that  obtained  by  Crawford  and 
Lemlich  (1962).  The  expansions  for  temperature  and  stream  function 
were  carried  out  to  three  terms.  It  was  estimated  that  for  Pr  ~  1 , 
the  convection  terms  were  negligible  in  comparison  to  the  conduction 
terms  for  Rayleigh  numbers  ranging  up  to  approximately  10  and  R  (the 
radius  ratio)  in  the  range  of  1.15  to  4.15.  Vertical  symmetry  was 
assumed,  and  for  Pr  >  1 ,  the  flow  resulted  in  the  single-cell  kidney¬ 
shaped  pattern  for  moderate  Rayleigh  numbers.  But  for  A  ~  0  (104),  the 


second  and  third  terms  in  the  stream  function  expansion  started  to 
outgrow  the  first  term,  so  any  resulting  flow  patterns  were  deemed 
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invalid.  For  the  case  of  Pr  =  .02  (liquid  mercury),  A  =  300  and 
R  =  2.0,  a  multicellular  flow  was  reported.  Two  weak  secondary  cells 
formed  at  the  top  and  bottom,  while  a  stronger  primary  cell  formed  near 
the  center  of  the  annulus.  At  the  time  though,  experimental  results  for 
the  larger  gap  low-Prandtl  number  range  were  unavailable;  hence,  their 
multicellular  flow  pattern  could  not  be  fully  supported.  Note  that  a 
similar  perturbation  analysis  was  performed  for  spherical  annuli  in 
Mack  and  Hardee  (1968).  Their  flow  field  patterns  were  very  similar 
to  those  obtained  in  the  concentric  cylinder  geometry. 

Rotem  (1972)  studied  the  conjugate  problem  of  conduction  within 
the  inner  cylinder  coupled  with  convective  motion  in  the  gap.  Through 
a  trial-and-error  procedure,  he  was  able  to  obtain  the  following 
expansions  for  stream  function  and  temperature: 

*  =  G  ^  +  G2^  +  Pr  G2^  +  G3^  +  Pr  G3^ 

+  Pr2  G3^  +  ....  (1  <  r  <  R)  (2.2a) 

T  =  To  +  RaT1  +  RAGT®  +  R2t£  +  RaGT®  +  ....  (2.2b) 

where  G  was  the  Grashof  number  based  on  (R-l),  RA  represented  the 
Rayleigh  number,  and  R  the  outer-to-inner  radius  ratio.  The  first  term 
in  these  expansions  corresponded  to  the  creeping-flow  solution. 

Rotem  (1972)  carried  out  the  stream  function  expansion  to  two  terms 
and  calculated  three  terms  for  the  temperature  perturbation  expansion. 


Through  his  results,  Rotem  (1972)  confirmed  the  basic  single-cell 
solution  and  reinforced  the  idea  that  one  should  be  careful  in  reporting 
counter-rotating  cells  with  an  asymptotic  method,  since  this  usually 
signifies  that  the  expansions  are  starting  to  diverge.  He  did  not 
investigate  the  extreme  Prandtl  number  cases,  but  he  did  express  the 
fact  that  further  transformations  were  needed  to  eliminate  Pr  as  an 
independent  parameter  and  render  the  equations  free  from  singularities 
in  the  limits  of  Pr  -*■  0  and  Pr  -+  ».  Such  has  been  obtained,  for 
example,  in  his  analysis  on  natural  convection  above  unconfined 
horizontal  surfaces  (Rotem  and  Claassen,  1969). 

Hodnett  (1973)  used  a  perturbation  method  to  analyze  the  same 
problem  as  Mack  and  Bishop  (1968),  except  that  his  analysis  was  in 
terms  of  primitive  variables.  He  extended  his  work  in  order  to  determine 
how  large  R  could  be,  at  a  given  value  of  Grashof  number,  for  the 
problem  to  remain  conduction  dominated.  He  found  that  convection  was 
negligible  only  when  R  satisified 

R3  Jin-1  R  =  0[(eG)_1]  (2.3) 

where  eG  represented  his  Grashof  number  and  e  was  given  by 


Huetz  and  Petit  (1974)  performed  a  theoretical  study  of  free 
convection  in  a  horizontal  annulus  for  low  values  of  Grashof  number. 
Their  governing  equations  were  written  in  the  same  form  as  Mack  and 
Bishop  (1968),  and  vertical  symmetry  was  also  assumed.  They  expanded 
stream  function  and  temperature  in  a  power  series  with  respect  to 
Grashof  number  (Gr),  where  Gr  was  based  on  the  inner  cylinder  radius. 

Two  case  studies  were  investigated.  Case  I  related  to  a  constant  heat 
flux  imposed  on  the  inner  wall  and  constant  temperature  on  the  outer 
wall,  and  vice-versa  for  case  II.  For  Pr  ~  1 ,  only  monocellular  flow 
was  obtained  in  both  cases  I  and  II,  regardless  of  the  radius  ratio  for 
Grashof  numbers  less  than  1,000.  However,  for  Pr  =  .02,  R  =  2.0,  and 
Gr  =  15,000  (or  R a  =  300),  a  multicellular  flow  was  observed.  Secondary 
cells  formed  at  the  top  and  bottom  of  the  annulus  with  the  primary  cell 
in  the  center,  but  the  secondary  cells  did  not  appear  simultaneously. 

In  case  II,  the  bottom  cell  appeared  first,  and  in  case  I,  the  upper 
cell  appeared  first.  Thus,  the  results  of  this  study  further  support 
the  multicellular  flow  recognized  by  Mack  and  Bishop  (1968). 

Custer  and  Shaughnessy  (1977a)  investigated  natural  convection 
within  a  horizontal  annuli  for  very  low  Prandtl  numbers  by  solving  the 
dimensionless  thermal  energy  and  vorticity  equations  with  a  double 
perturbation  expansion  in  powers  of  Grashof  and  Prandtl  numbers.  The 
stream  function  and  temperature  expansions  were  written  as: 


(2.4b 


T(r,8)  «  E  E  PrJ  Gr  T..  (r,0) 
j=0  k=0  o  JK 

For  a  fairly  large  gap  of  R  =  5  and  Pr  -*■  0,  they  reported  that  the 

center  eddy  fell  downward  as  the  Grashof  number  increased,  which  was 

contrary  to  the  behavior  of  fluids  for  Pr  >_  .7.  For  this  same  size 

gap,  at  Gr  =  12,000,  they  observed  the  formation  of  a  weak  eddy  near 
ro 

the  top  of  the  inner  cylinder.  For  R  =  2,  at  GrrQ  =  120,000,  two  weak 
eddies  formed  near  the  top  and  bottom  of  the  annulus  while  the 
stronger  kidney-shaped  one  remained  in  the  center.  Vertical  symmetry 
was  again  assumed  for  this  problem.  These  conclusions  also  help  to 
confirm  the  multicellular  flow  pattern  observed  by  Mack  and  Bishop 
(1968).  Custer  and  Shaughnessy  also  stressed  that  only  numerical 
solutions  to  the  full  nonlinear  equations,  or  experiments,  could 
actually  establish  the  true  existence  of  this  multicellular  flow  field 
Custer  and  Shaughnessy  (1977b)  continued  to  study  the  problem 
described  in  the  above  review.  In  this  analysis,  the  dependent 
variables  in  the  governing  equations  were  represented  by  the  following 
partial  spectral  expansions: 


oo 


^(r,y)  =  z  f  (r)  sin  n6 
i  n 
n=1 

(2.5a 

oc 

T(r,0)  =  2  g_(r)  cos  ne  . 

(2.5b 
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The  use  of  these  expansions  resulted  in  ordinary  differential  equations 
governing  f  (r)  and  gn(r).  The  equations  were  solved  numerically  and 
their  results  were  based  on  the  series  being  truncated  after  three 
terms.  For  Pr  =  .01  and  Ra  <  50  (or  Gr  <  5,000),  they  obtained  the 
following  correlation  for  their  heat  transfer  data  (note  that  Ra  and  Gr 
were  based  on  the  inner  cylinder  radius): 


at  R  =  2.00, 

keq  =  1  +  5  x  10"8  Ra2,014 

(2.6a) 

at  R  =  4.00, 

keq  =  1  +  1  x  10"4  Ra1,935  , 

(2.6b) 

where  k&q  signifies  the  average  equivalent  thermal  conductivity  and 
represents  the  actual  mean  heat  flux  divided  by  the  heat  flux  for  pure 
conduction.  For  a  radius  ratio  of  1.1,  they  found  that  kgq  was  equal 
to  1  out  to  a  Grashof  number  of  approximately  106.  Also,  for  all  the 
cases  studied,  only  the  single-cell  flow  field  resulted  in  the  vertical 
half  of  the  annulus. 

Walton  (1980)  utilized  a  multiple-scales  linearized  stability 
theory  to  study  the  instability  of  natural  convective  flow  between 
narrow  cylindrical  annuli.  First,  he  represented  the  basic  flow  by 
expanding  the  stream  function  (>')  and  temperature  (T)  in  terms  of  the 
dimensionless  gap  width,  e: 


where 

r  -  r. 

e  =  — -r~- 1  and  e  «  1  .  (2.7b) 

He  then  let  and  T'  represent  small  perturbations  to  and  T  given 

in  Eq.  (2.7a).  By  using  the  dimensionless  forms  of  the  vorticity  and 
energy  equations,  he  was  able  to  consider  the  stability  of  the  basic 
flow  to  small  disturbances,  via  linearized  stability  theory.  He 
determined  that  the  convective  flow  became  unstable  at  a  critical  value 
of  the  Rayleigh  number,  R,  (associated  with  a  particular  wavenumber) 
given  by: 

R  =  Rc  +  e  R1  +  ....  (2.8) 

where,  for  Pr  =  .7,  Rc  =  1707.762  and  R^  =  258.4.  Thus,  as  e  0, 

R  1707.762.  According  to  his  results,  the  narrow-gap  annulus 
collapsed  to  the  horizontal  flat  plate  Benard  problem  as  the  gap  width 
tended  to  zero.  Hence,  for  Pr  ~  1 ,  a  thermal -type  instability  should 
arise  near  the  top  of  the  annulus. 

Jischke  and  Farshchi  (1980)  studied  the  boundary-layer  regime  for 
laminar  free  convection  between  2-D  horizontal  annuli  at  the  large 
Rayleigh  number  limit.  They  divided  the  flow  field  into  five 
physically  distinct  regions,  valid  for  the  high  Rayleigh  number 
limiting  condition.  They  assumed  a  stagnant  regime  for  the  bottom  part 
of  the  annulus;  a  boundary-layer  type  behavior  near  the  inner  and 
outer  cylinders;  an  inviscid  core  region  in  the  center  portion  of  the 


12 


annulus;  plus  a  thermal  plume  section  along  the  vertical  line  of 
symmetry  above  the  inner  cylinder,  where  the  inner  boundary-layer  joins 
up  with  the  outer  boundary-layer.  Their  governing  equations  were 
written  in  terms  of  primitive  variables,  and  they  employed  a  zeroth-order 
asymptotic  expansion  to  represent  the  velocity  and  temperature  fields 
within  the  annulus.  For  the  boundary-layer  regime,  they  assumed  (from 
free  convective  boundary-layer  studies  of  a  single  horizontal  cylinder) 
that  the  velocity  and  temperature  fields  would  scale  as: 

u  ~  Ra1/4  uQ  +  _ 

v  ~  v  +  . . . . 
o 

T  ~  T0  +  ....  (2.9) 

where  u  signified  the  radial  velocity  and,  v,  the  tangential  velocity 
component.  Their  simplified  equations,  resulting  from  the  Ra  ® 
limit,  were  solved  by  means  of  an  integral  method  in  the  limit  of 
Pr  °o.  They  compared  their  heat  transfer  results  to  those  of  Kuehn 
and  Goldstein  (1976a)  for  a  radius  ratio  of  R  -  2.6,  Pr  =  .706  and 

4 

Ra  =  4.7  x  10  .  Qualitatively,  their  basic  flow  field  and  results  were 
very  similar,  but  significant  deviation  was  apparent  in  their  plots  of 
the  variation  of  local  Nusselt  number  with  angular  position  on  the 
inner  cylinder.  Their  best  agreement  was  achieved  near  the  top  of  the 
inner  cylinder,  at  9  =  180°.  They  seemed  to  have  captured  the 
essential  features  of  the  flow  field  in  their  boundary-layer  analysis, 


although  only  bicel lular-type  flows  were  studied.  The  authors  of  this 
article  did  not  investigate  the  narrow-gap  limit  of  their  Ra  -*■  °° 
equations. 

2.2.  Experimental  Studies 

Most  of  the  articles  in  this  section  relate  to  natural  convection 
between  either  concentric  cylinders  or  spheres.  Various  flow  patterns 
have  been  investigated  in  the  experimental  work,  ranging  from  the 
multicellular  flow  prevalent  in  the  narrow  gaps  to  the  unsteady  type 
flows  originating  in  the  larger  gaps.  Studies  include  both  2-D  and  3-D 
phenomena  over  a  range  of  fluids  such  as  air,  water,  liquid  mercury,  oil 
and  glycerin. 

Liu  et  al_.  (1962)  studied  natural  convection  heat  transfer  for  air, 
water  and  silicone  oil  (0.7  ±  Pr  3,500)  between  horizontal 
cylindrical  annuli.  They  considered  five  different  geometries  with 
radius  ratios  of  R  =  1.154,  1.5,  2.5,  3.75  and  7.5.  For  the  larger 
radius  ratios,  all  three  fluids  experienced  a  slow  sideways  oscillation 
near  the  top  of  the  cylinders  as  the  Rayleigh  number  (Ra.  )  exceeded 

u-a 

some  limiting  value,  and  this  valt.-.  decreased  with  gap  size.  For  the 
smaller  gap,  R  =  1.154,  a  multicellular  type  flow  (near  the  top)  was 
observed  for  air  and  silicone  oil,  at  Ra^_a  of  approximately  2,000 
and  18,000,  respectively.  A  kidney-shaped  flow  pattern  was  maintained 
below  the  counter-rotating  cells  that  formed  near  the  top  of  the 
annulus.  They  speculated  that  for  gaps  smaller  than  R  =  1.154,  the 
Benard  critical  Rayleigh  number  (Ra.  *  1700)  would  be  approached. 


Bishop  ejt  aj_.  (1964)  investigated  natural  convective  heat  transfer 
between  concentric  spheres.  They  examined  two  radius  ratios,  R  =  3.14 
and  R  =  1.19.  For  the  larger  gap,  only  the  kidney-shaped  flow  pattern, 
similar  to  that  seen  for  concentric  cylinders,  was  observed.  No 
sideways  oscillations  were  noticed.  But  in  the  smaller  gap,  R  =  1.19, 
two  counter-rotating  cells  appeared  near  the  top  at  a  Rab  of  about 
3,600.  However,  due  to  the  spherical  geometry,  the  two  secondary  cells 
almost  immediately  began  to  coalesce  into  an  elongated  shape  and  soon 
became  indistinguishable.  Although  the  geometries  are  different, 
the  multicellular  flow  observed  in  the  narrow  spherical  annuli  supports 
the  same  type  of  flow  seen  by  Liu  et^  al_.  (1962)  in  the  narrow 
cylindrical  annuli. 

Grigull  and  Hauf  (1966)  used  a  Mach-Zehnder  interferometer  to 
measure  the  temperature  field  between  horizontal  cylindrical  annuli 
filled  with  air.  From  these  measurements,  they  were  able  to  obtain  the 
local  Nusselt  numbers  for  a  particular  gap  size  and  Grashof  number.  They 
discussed  three  different  regimes  of  convective  flow: 

1.  A  2-D  pseudo-conductive  regime  for  Gr.  <  2,400.  Here, 
conduction  effects  were  dominant,  although  some  convective 
motion  was  evident. 

2.  A  transitional  regime  with  3-D  convective  motion,  for 
2,400  <  Gr,  ,  <  30,000.  For  the  intermediate-size  gaps, 

(1.2  <  R  <  2.0),  the  flow  field  transitioned  to  a  form  of 
3-0  vortices  in  unsteady  oscillatory  motion. 


M-" "■’’v »«■  V v ^ v." TOT Vffl" TO BWV ITMW  V  VAT  v  ^A^vv^.wmw.WTi 


15 

3.  A  2-D  fully-developed  laminar  convective  motion  regime 

for  30,000  <  Gr,  _  <  716,000.  The  2-D  motion  was  considered 

—  D  *  a  — 

steady  for  the  larger  gap  widths  in  this  particular  range  of 
Grashof  numbers.  They  did  notice  that  as  Grashof  number 
increased,  the  centers  of  the  kidney-shaped  cells  moved 
closer  to  the  upper  portion  of  the  annulus. 

Vivid  flow  pictures  (observed  by  cigarette  smoke)  were  included  for  all 
three  convective  flow  regimes. 

Bishop  et  al_.  (1966)  extended  their  study  of  1964  (Bishop  et  al . , 
1964)  to  include  the  effects  of  various  radius  ratios,  ranging  from 
R  =  1.19  to  R  =  3.14.  In  their  experiments,  they  confirmed  the 
occurrence  of  counter-rotating  cells  previously  observed  in  the  smallest 
gap  width,  R  =  1.19.  For  the  relatively  narrow  gaps  (slightly  greater 
than  R  =  1.19)  in  the  high  Grashof  number  range,  they  observed  a 
boundary-layer  type  flow  near  the  walls  of  the  annulus,  coupled  with  a 
slower,  more  uniform  type  of  fluid  motion  near  the  center.  Also,  from 
their  measured  heat-transfer  data  for  four  different  radius  ratios, 
they  obtained  two  Nussel t-Grashof  number  correlations  that  fit  their 
data  to  within  15.5  percent. 

Lis  (1966)  studied  the  flow  behavior  in  simple  and  obstructed 
annuli  using  the  Schileren  technique  with  both  sulphur  hexafluoride 
and  nitrogen  as  the  working  fluid.  The  six  axial  spacers  used  in  the 
obstructed  annuli  provided  for  enhanced  heat  transfer  effects  (compared 
to  the  simple  annuli)  due  to  more  efficient  mixing  near  the  upper  parts 
of  the  annulus.  Also,  for  2  <  R  <  4,  the  flow  pattern  became  unstable 


at  the  higher  Rayleigh  numbers  and  random  oscillations  developed  near 
the  top  of  the  gap  (for  the  unobstructed  annuli). 

Bishop  and  Carley  (1966)  performed  photographic  studies  of  natural 
convective  air  patterns  between  concentric  cylinders.  They  observed  an 
oscillatory  flow  for  R  =  3.69,  that  started  at  Ra.  =  270,000  (a  much 
larger  value  than  that  indicated  by  Liu  et  al_.  (1962)  for  the  same  size 
gap).  They  also  tried  to  reproduce  the  multicellular  flow  field 
observed  by  Liu  et  al_.  (1962)  for  the  small  gap,  R  =  1.154.  They  used 
a  gap  width  of  .688  inches,  and  for  Rayleigh  numbers  (Ra,  )  up  to 
20,000,  no  type  of  multicellular  flow  occurred.  Either  the  gap  size 
was  too  large,  or  a  much  higher  Rayleigh  number  was  needed  to  observe 
the  cells  with  their  particular  apparatus. 

Bishop  et  al_.  (1968)  re-examined  experimentally,  for  air,  the 
oscillatory  flow  in  the  larger  size  cylindrical  annuli.  For  gap  widths 
of  R  =  3.70,  2.46  and  1.846,  they  observed  that  the  oscillations 
started  at  Rayleigh  numbers  (Ra,  _)  of  approximately  240,000,  150,000 
and  35,000,  respectively.  Using  the  data  from  these  three  gaps,  they 
obtained  empirical  correlations  for  the  period,  wavelength  and 
amplitude  of  the  oscillations. 

Powe  et  a]_.  (1969)  experimentally  studied  the  natural  convective 
flow  of  air  between  horizontal  concentric  cylinders  for  various  radius 
ratios.  They  characterized  the  flow  into  three  basic  regimes: 

I.  2-D  multicellular  flow  for  R  <  1.2, 

II.  3-D  spiral  flow  for  1.2  <  R  <_  1.7,  and 

III.  2-0/3-0  oscillatory  flow  for  R  >  1.7. 
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Regimes  II  and  III  were  unsteady.  The  multicellular  flow  field  in 
Regime  I  remained  steady  for  a  small  number  of  cells,  but  as  the 
Rayleigh  number  increased  and  the  chain  of  cells  near  the  top  stopped 
forming,  slight  oscillations  about  the  vertical  center-line  occurred. 

Yin  et  al .  (1973)  experimental ly  investigated  the  natural 
convective  flow  patterns  between  isothermal  concentric  spheres  for  air 
and  water.  For  air,  with  R  =  1.4  and  Rab_a  =  5930,  they  observed  two 
steady  counter-rotating  cells  near  the  top  of  the  annulus.  With  the 
larger  gaps  (R  =  2.17  and  1.78),  for  both  air  and  water,  the  familiar 
kidney-shaped  pattern  resulted  for  relatively  small  Rayleigh  numbers. 
But  as  Rayleigh  number  increased,  an  unsteady  flow  behavior  was 
initiated. 

Kuehn  and  Goldstein  (1976a)  studied  the  natural  convective  flow  of 
air  and  water  between  isothermal  concentric  cylinders  both 
experimentally  and  numerically.  A  Mach-Zehnder  interferometer  was  used 
to  obtain  temperature  profiles  for  both  air  and  water  with  an  annulus 
of  R  =  2.6.  Heat  transfer  correlations  were  found  by  using  a  least- 
squares  regression  analysis.  For  air, 

k  =  .159  Ra’2^2  ,  2.1  x  104  <  Ra.  <  9.6  x  104  (2.10) 

eq  b-a  b-a 

and  for  water, 

k  =  .234  Ra/288  ,  2.3  x  104  <  Ra.  <  9.8  x  105  .  (2.11) 

eq  b-a  b-a 


At  similar  Rayleigh  numbers,  the  temperature  distributions  for  air  and 
water  were  essentially  the  same.  The  flow  field  remained  steady  and 
symmetric  for  all  Rayleigh  numbers  investigated  at  R  =  2.6.  There  was 
no  sign  of  any  type  of  multicellular  flow  behavior. 

Kuehn  and  Goldstein  (1978)  studied  the  effects  of  eccentricity  and 

Rayleigh  number  on  natural  convection  heat  transfer  between  horizontal 

annuli  filled  with  pressurized  nitrogen.  A  Mach-Zehnder  interferometer 

was  again  used  to  obtain  interferograms  of  the  working  fluid.  For  an 

2 

eccentric  geometry  of  ev/L  <_^  (meaning  the  inner  cylinder  is  set 
above  the  concentric  center),  the  overall  heat  transfer  rates  were 
within  ten  percent  of  that  for  concentric  cylinders  at  the  same  Rayleigh 
number.  But  local  changes  in  the  heat  transfer  rates  were  significant 
where  the  cylinder  walls  were  nearly  touching.  When  the  inner  cylinder 
was  closest  to  the  bottom  of  the  outer  cylinder  (ev/L  -  -.623),  an 
unsteady  thermal  plume  behavior  was  witnessed,  with  the  flow  eventually 
transitioning  to  turbulence  upon  further  increase  in  Rayleigh  number. 

Warrington  and  Powe  (1985)  experimentally  examined  natural 
convection  between  concentrically  located  isothermal  spherical, 
cylindrical  and  cubical  inner  bodies,  surrounded  by  a  cubical 
enclosure.  To  within  15  percent  of  their  experimental  data,  mean 
Nusselt  number  correlations  were  obtained  for  all  three 
configurations.  In  comparison  to  heat  transfer  data  from  spherical 
annuli,  they  found  that  cubical  enclosures  with  a  spherical  inner 
body  yielded  larger  Nusselt  numbers  for  a  given  Rayleigh  number 
and  aspect  ratio.  Also,  for  the  higher  Rayleigh  numbers,  a 


multi  cellular  flow  field  was  encountered  for  the  cubical  enclosures 


with  cylindrical  inner  bodies. 

2.2.1.  Heat  transfer  correlations 

Some  basic  correlations  are  presented  in  McAdams  (1954)  for 
natural  convection  from  single  horizontal  cylinders.  He  observed  that 
the  heat  transfer  data  for  single  cylinders  should  approach  that 
obtained  with  concentric  cylinders  in  the  limit  as  R  •+■  °°.  McAdams 
also  presented  expressions  for  natural  convective  heat  transfer  between 
horizontal  and  vertical  enclosed  air  spaces.  All  of  his  correlations 
took  on  the  form  Nu  =  a(Ra)b,  where  b  was  typically  1/4  for  laminar 
flows  and  1/3  for  turbulent- type  flows. 

Itoh  et  aJL  (1970)  proposed  a  new  method  for  correlating  heat 
transfer  coefficients  for  natural  convection  between  horizontal 
cylindrical  annuli.  They  claimed  that  heat  transfer  coefficients  (hi , 
hQ)  are  well  correlated  by  the  mean  Nusselt  number,  NlT,  and  the  mean 
Grashof  number,  Gr^,  defined  as  follows: 

^  h  -Cr  S.n(r  /r  )]  h  -[r  ln(r  /r  )] 

Nu  ■ - k - =  - k -  <2J2> 

96(1,  -  T,)[(/r,r  )nn(i-  /r,)]3 

Grm  =  - i - 2 - J-2 - - — - -  .  (2.13) 

V 

With  these  definitions,  any  resulting  heat  transfer  data  (for  laminar 
flow)  could  be  represented  by 

"\  /  A 
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where  c-|  is  a  constant  value. 

Powe  (1974)  examined  heat  transfer  correlations  given  by  Liu  et  al . 
(1962)  and  Scanlan  et^  al_.  (1970)  to  determine  the  bounding  effects  of 
heat  loss  by  free  convection  from  concentric  cylinders  and  spheres, 
respectively.  For  various  Rayleigh  numbers,  he  calculated  limiting 
bands  at  which  the  empirical  equations  were  valid.  Beyond  these  limits, 
the  equations  either  collapsed  to  the  conduction  solution  for  small  gap 
widths,  or  to  natural  convective  flow  of  a  single  cylinder  (or  sphere) 
exposed  to  an  infinite  atmosphere,  for  the  larger  gap  sizes. 

Raithby  and  Hollands  (1975)  proposed  a  heat  transfer  correlation 
for  isothermal  concentric  cylinders.  Based  on  their  experimental  data 
for  air,  water  and  silicon  oil,  they  obtained  the  following  correlation 
(valid  for  the  convection-dominated  flow  regime): 

keq  =  .386  [Pr/( .861  +  Pr)]1/4  Rac(.1/4  (2.14b) 

where 


O  (D  /D  )]4 

Ra  =  - 5 - X-TF - 0  ,-r  -r  Ra.  , 

cc  (b-a)3  (1/Di 3/5  +  1/Dq3/5)5  b 


(2.14c) 


and  b-a  =  annulus  gap-width.  They  mentioned  that  the  correlation 
worsened  as  the  annular  gap  spacing  increased,  but  appeared  to  be 
highly  satisfactory  for  the  relatively  narrow  and  intermediate  size 
gaps. 
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Kuehn  and  Goldstein  (1976b)  obtained  correlating  equations  for 
heat  transfer  between  horizontal  concentric  cylinders.  For  R  =  2.6 
and  .01  £  Pr  <  1,000,  the  following  expression  (for  the  mean  inner- 
diameter  Nusselt  number)  represented  their  heat  transfer  data  to 
within  2  percent: 


Nu 


D,- 


in 


(([( — "~25-)5/3+(-587  G  Ra|!/4)5/3]3/5)l5+( .  1  Ra^3)15}1715 


1-e 


(2.15) 

where 

G  =  [(1  +  9-^L)  +  (.4  +  2.6  Pr,7)"5]'1/5 

Pr’7  I 

I 

and  Ra^  ,  RaQ  signify  the  Rayleigh  numbers  evaluated  at  the  inner  ! 

i  o  I 

i 

and  outer  cylinder  diameters.  They  also  discussed  that  as  the  diameter  < 

t 

of  the  outer  cylinder  increased,  the  heat  transfer  approached  that  of  j 

i 

a  single  horizontal  cylinder.  They  found  that  to  have  heat  transfer 
within  5  percent  of  a  free  cylinder  required  DQ/ >  360  at 
RaQ  =  107  and  Dq/D.  >  700  for  Ra^  =  10'1 
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Boyd  (1981)  used  a  unified  theory  to  correlate  steady  laminar 


natural  convection  heat  transfer  between  horizontal  annuli.  His 
correlations  were  successfully  extended  to  annuli  with  irregular 
boundaries.  For  isothermal  concentric  cylinders,  he  suggested  the 
following  relationship  for  the  mean  Nusselt  number  based  on  the  gap 
width.  A: 


Nua  «  C4  Prn*  Ra1/4  (2.16) 

where 

n*  =  C5  +  C6  Pr'1/3 

and  C^,  and  were  constants  that  depended  on  related  heat  transfer 
data.  These  constants  were  evaluated  using  data  from  Keuhn  and 
Goldstein  (1976a).  The  result  is  given  by 

Nu.  =  .796  Prn*  Ra1/4  (2.17) 

A 


where 

n*  =  .00663  -  .0351  Pr'1/3  . 

According  to  Boyd,  this  expression  was  valid  for  10^  <  Ra  <  107, 

.706  1  Pr  3100,  and  .125  <_  —  <_  2.0;  where  represented  the  aspect 

ri  ri 


ratio  for  the  annulus. 
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2.3.  Numerical  Studies 

As  previously  mentioned,  most  of  the  numerical  studies  for 
natural  convection  between  horizontal  annuli  pertain  to  the  bicellular 
solutions  associated  with  the  low  to  moderate  Rayleigh  numbers. 

Although  the  reviews  in  this  section  will  be  rather  concise,  particular 
attention  will  focus  on  the  finite-differencing  methods  used,  especially 
with  regard  to  the  numerical  representation  of  the  nonlinear  convective 
terms. 

Crawford  and  Lemlich  (1962)  studied  natural  convection  of  air 
between  horizontal  cylindrical  annuli.  They  numerically  examined  three 
different  radius  ratios,  R  =  2,  8,  and  57,  and  confined  their  study 
to  extremely  low  Grashof  numbers,  the  so-called  creeping  flow  solution. 
In  their  numerical  method,  conventional  central -di fferencing  was  used 
throughout,  and  the  stream  function  and  temperature  were  calculated 
using  a  Gauss-Seidel  iterative  procedure.  Vertical  symmetry  was 
assumed  and  their  results  revealed  the  characteristic  kidney-shaped 
circulation  pattern. 

Abbott  (1964)  discussed  a  numerical  method  for  solving  the  same 
problem  as  Crawford  and  Lemlich  (1962),  except  for  very  narrow  annuli. 

He  studied  four  types  of  narrow  gaps,  R  =  1.0256,  1.0170,  1.0084  and 
1.0040.  Abbott  began  his  solution  process  by  first  obtaining  solutions 
to  the  conduction-dominated  energy  equation  and  the  creeping-flow 
(negligible  convective  acceleration  terms)  momentum  equation.  He  then 
used  these  results  to  approximate  the  convective  terms  in  the  full  set 
of  governing  equations.  Thus,  by  linearizing  the  equations  in  terms 
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of  the  unknowns,  he  was  able  to  obtain  successive  approximations  to  the 
nonlinear  terms.  He  solved  for  the  unknown  temperatures  and  stream 
functions  by  using  a  matrix  inversion  technique.  His  results  represent 
only  slight  convective  perturbations  to  the  creeping-flow  solution. 

For  all  cases,  he  obtained  the  basic  kidney-shaped  flow  pattern. 

Powe  et^  al_.  (1971)  obtained  the  first  semblance  of  a  secondary 
air  flow  captured  numerically  in  narrow  horizontal  annuli.  They  assumed 
vertical  syrnnetry  and  used  a  complete  central -difference  representation 
for  the  nonlinear  convective  terms  in  the  governing  energy  and 
vortici ty-transport  equations.  Their  equations  were  formulated  for 
the  steady-state  case.  Moreover,  their  numerical  solutions  could  only 
capture  the  flow  field  up  to  the  point  where  the  stream  function  changed 
its  sign  (signifying  counter-rotating  flow).  After  this  point,  their 
code  would  no  longer  converge.  Thus,  they  could  not  fully  resolve  the 
mul ticel 1 ular  flow  field  for  the  narrow  gaps,  but  they  were  able  to 
make  estimates  of  the  transitional  Rayleigh  numbers.  Using  this 
approach,  they  estimated  a  transitional  Rayleigh  number  of 
Ra  =  452,000  for  R  =  1.2.  Their  mean  Nusselt  numbers  were  not  at  all 
affected  by  the  appearance  of  these  unresolved  secondary  flows. 

Charrier-Motjtabi  et  al_.  (1979)  used  an  ADI  scheme  to  numerically 
solve  for  the  natural  convective  flow  field  between  horizontal  annuli. 
Coupled  with  the  energy  equation,  they  used  the  stream  function- 
vorticity  approach  to  analyze  the  flow.  Also,  a  fictitious  time  was 
defined  so  as  to  more  readily  achieve  the  steady-state  condition.  They 
did  not  mention  what  type  of  finite-differencing  was  employed  for  the 
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nonlinear  convective  terms.  For  the  narrow  gap  of  R  =  1.2  and  Rayleigh 
numbers  up  to  Ra  =  875,000,  no  multicellular  flow  field  was  observed 
for  air.  However,  for  Pr  =  .02  and  R  =  2.0,  they  obtained  a 
multicellular  flow  similar  to  that  described  by  Mack  and  Bishop  (1968). 

Asti  1 1  et  al_.  (1979)  obtained  numerical  solutions  for  natural 
convection  in  concentric  spherical  annuli.  They  considered  fluids  with 
Prandtl  numbers  between  .7  and  5.00,  and  radius  ratios  from  1.03  to 
2.00.  Pure  central -differencing  was  used  to  obtain  approximations  to 
the  steady-state  stream  function  and  energy  equations.  The  system  of 
equations  was  solved  using  a  Gauss-Seidel  iterative  procedure  with 
under-relaxation,  and  vertical  symmetry  was  again  assumed.  They  mainly 
observed  the  typical  kidney-shaped  pattern  for  all  fluids  and  gaps 
studied.  However,  for  air  and  R  =  1.1  and  1.2,  they  resolved  a  second 
vortex  formation  occurring  at  Ra^  =  25,000  and  8,000,  respecti vely . 
These  counter-rotating  cells  were  seen  experimentally  by  Bishop  et  al . 
(1964),  but  at  a  much  lower  Rayleigh  number  for  R  =  1.2. 

Caltagirone  et  al_.  (1979)  also  performed  a  numerical  study  of 
natural  convection  between  spherical  annuli  filled  with  air.  They 
used  an  ADI  scheme  to  solve  for  the  stream  function,  temperature  and 
vorticity  fields.  They  considered  vertical  symmetry  and  radius  ratios 
between  1.15  and  3.00.  Using  "zero"  initial  conditions,  only 

5 

bicellular  solutions  were  found  for  Rayleigh  numbers  up  to  10  .  But 
upon  using  an  initial  temperature  distribution  of  the  form 


T  <*  a  sin  (-nr)  cos  (be) 


(2.18) 
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(where  'a*  was  an  amplitude  value  and  'b1  a  characteristic  wavenumber), 
the  converged  flow  field  for  R  =  2  and  Ra  =  50,000  experienced  a 
counter-rotating  cell  near  the  top  of  the  half-annulus.  Associated 
with  this  secondary  flow  was  a  rise  in  the  mean  Nusselt  number. 

But,  this  multicellular  flow  field  was  probably  unrealistic,  since 
this  type  of  flow  has  only  been  observed  experimentally  for  the  very 
narrow  spherical  annuli,  R  £  1.4. 

Kuehn  and  Goldstein  (1980b)  studied  the  effects  of  Prandtl  number 
and  diameter  ratio  on  natural  convection  between  horizontal  cylindrical 
annuli.  They  employed  an  explicit  finite-difference  scheme  to  solve 
the  energy,  stream  function  and  vorticity  equations  for  steady  laminar 
flow.  They  used  a  hybrid  technique  for  the  nonlinear  terms,  which 
switched  from  central  to  upwind  differencing  when  the  mesh  Reynolds 
number  constraint  was  exceeded.  For  large  Prandtl  numbers,  a  fully- 
developed  boundary-layer  with  an  impinging  thermal  plume  resulted, 
whereas,  as  Pr  -*•  0,  the  temperature  distribution  approached  the  pure 
conduction  limit.  Also,  the  mean  Nusselt  number  asymptotical ly 
approached  the  single  horizontal  cylinder  value  as  R  ->•  °°.  The  authors 
derived  a  correlation  valid  for  laminar  flow  over  .001  £  Pr  £  1,000 
and  1-0  £R  <  00  (although  the  correlation  fit  best  for  Pr  =  .7): 
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where  RaQ  and  RaQ  are  the  same  as  those  described  in  Eq.  (2.15). 
ui  uo 

Over  their  range  of  data,  they  reported  a  maximum  deviation  of  13 
percent  occurring  for  Prandtl  numbers  near  .02. 

Projahn  et  aj_.  (1981)  solved  the  energy  and  the  vortici ty-stream 
function  equations  numerically  by  using  a  strongly  implicit  method  as 
described  by  Weinstein  e_t  aj_.  (1970).  The  convective  terms  were  written 
in  divergence  form  and  were  differenced  using  a  corrected  upwind  scheme 
obtained  from  Jacobs  (1973).  In  their  analysis  of  natural  convection 
between  concentric  and  eccentric  cylinders,  they  reported  that  for  a 
negative  (downward)  vertical  displacement  of  the  inner  cylinder,  the 
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mean  Nusselt  number  was  always  greater  than  that  in  the  concentric  case. 
Also,  for  Pr  =  .7  and  R  =  2.6,  they  obtained  a  counter-rotating  cell 
near  Ra  =  12,000,  but  attributed  this  cell  formation  to  their 
assumption  of  symmetrical  boundary  conditions. 

Ingham  (1981)  solved  a  set  of  equations  similar  to  that  described 
by  Projahn  et_  al_.  (1981),  but  Ingham's  equations  were  formulated  in  a 
more  general  fashion  to  account  for  either  the  concentric  spherical  or 
cylindrical  geometry.  He  assumed  steady-state  and  vertical  symmetry 
and  used  central -differencing  throughout,  except  for  the  nonlinear 
terms  where  he  employed  a  cleverly  weighted  second-order  upwind- 
differencing  scheme.  He  tested  for  steady  multicellular  flows  by 
considering  a  spherical  radius  ratio  of  R  =  1.19  with  Pr  =  .7,  and 
for  a  range  of  values  of  Ra  up  to  2.5  x  107,  he  did  not  obtain  any  sign 
of  a  multicellular  structure. 

Farouk  and  Guceri  (1982)  were  the  first  to  study  turbulent  natural 
convection  numerically  between  2-D  horizontal  concentric  cylinders. 

They  used  a  k  -  e  turbulence  model  and  obtained  steady-state  results 
that  were  in  good  agreement  with  experimental  data.  All  of  their 
test  cases  pertained  to  R  =  2.6  with  vertical  symmetry  assumed,  and 
they  considered  Rayleigh  numbers  (based  on  gap  width)  up  to  10^  -  107. 

Cho  et^  al_.  (1982)  studied  the  natural  convection  of  air  in 
eccentric  horizontal  isothermal  cylindrical  annuli.  They  solved  the 
problem  using  bipolar  coordinates  with  the  assumption  of  vertical 
symmetry.  Central -differencing  was  employed  for  all  the  derivatives, 
including  the  nonlinear  terms.  They  numerically  investigated  a 
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radius  ratio  of  R  =  2.6  with  varying  degrees  of  eccentricity,  for 
Raleigh  numbers  (based  on  gap  width)  less  than  5.0  x  10^.  Their 
numerical  results  agreed  rather  well  with  experimental ly  obtained 
interferograms.  It  was  observed  that  the  average  heat  transfer 
increased  as  the  inner  cylinder  moved  downward  along  the  vertical 
center-1 ine. 

Prusa  and  Yao  (1983)  numerically  examined  natural  convection 
between  eccentric  horizontal  cylinders,  similar  to  that  described  in 
the  above  with  Cho  et  al_.  (1982).  However,  they  employed  a  unique 
radial  transformation  method  that  allowed  them  to  study  various 
eccentricities  while  avoiding  any  type  of  singular  behavior  in  the 
limit  of  zero  eccentricity.  They  also  developed  a  convenient  variable 
mesh  routine  which  provided  the  flexibility  of  concentrating  grid 
nodes  within  the  boundary- layer  and  thermal  plume  regions.  Central- 
differencing  was  used  throughout,  along  with  a  stable  corrected  second- 
order  central  difference  scheme  to  represent  the  nonlinear  terms.  They 
obtained  very  good  agreement  with  experimental  and  analytical  data,  and 
confirmed  the  fact  that  the  overall  heat  transfer  could  be  reduced  or 
enhanced  with  respect  to  the  upward  or  downward  shift  of  the  inner 
cylinder  about  the  vertical  center-line.  In  addition,  they  appeared  to 
be  the  first  to  determine  the  critical  eccentricity  associated  with 
minimum  heat  transfer  for  various  Grashof  numbers. 

Chandrashekar  et_  al_.  (1984)  studied  natural  convective  flow  of  a 
Boussinesq  heat-generating  fluid  between  two  horizontal  concentric 
cylinders.  They  investigated  the  effects  of  two  driving  mechanisms  - 


an  externally  imposed  temperature  gradient  across  the  annulus,  coupled 
with  a  uniform  internal  heat  generation.  The  parameter  that 
represented  the  ratio  of  the  internal  heating  to  the  applied  temperature 
difference  was  denoted  by  S,  where  the  S  =  0  limit  corresponded  to 
the  isothermal  concentric  cylinder  case.  They  assumed  vertical 
symmetry  and  marched  numerically  in  time  to  achieve  the  steady-state 
condition.  The  governing  equations  were  solved  with  an  ADI  scheme  on 
a  uniform  mesh,  and  central -differencing  was  used  throughout.  They 
found  that  as  S  increased  from  zero,  a  transition  took  place  (for 
Pr  =  .7)  from  a  unicellular  to  a  bicellular  circulation  in  each  half¬ 
cavity.  For  this  transition,  the  critical  value  of  S  depended  on  both 
the  Rayleigh  number  and  the  radius  ratio. 

Tsui  and  Tremblay  (1984)  used  an  unsteady  code  to  obtain  steady- 
state  solutions  for  natural  convective  air  flow  between  horizontal 
annuli.  They  employed  an  ADI  scheme  to  solve  the  energy  and 
vorticity-stream  function  equations.  They  assumed  vertical  symmetry 
and  used  a  relatively  coarse  mesh  of  16  radial  nodes  together  with 
21  angular  nodes.  For  moderate  size  Rayleigh  numbers,  they  achieved 
steady-state  results  for  three  radius  ratios,  R  =  1.2,  1.5  and  2.0. 

In  all  of  these  cases,  they  obtained  the  steady-state  kidney-shaped 
flow  pattern. 

Lee  (1984)  studied  laminar  convection  of  air  between  concentric 
and  eccentric  heated  rotating  cylinders.  His  numerical  method 
involved  a  mesh  transformation  technique  coupled  with  the  introduction 
of  false  transient  time  terms  that  facilitated  steady-state  solutions 
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to  the  problem.  He  studied  a  radius  ratio  of  2.6  at  different  Rayleigh 
numbers,  Reynolds  numbers,  and  eccentricities.  His  governing  equations 
were  solved  with  an  ADI  method  and  all  spatial  derivatives  were 
approximated  by  second-order  central -differences.  Second-order 
upwind-differencing  was  used  on  the  convective  terms.  For  Ra  =  25,000, 
he  obtained  a  multicellular  type  flow  when  the  inner  cylinder  was 
shifted  upward  next  to  the  outer  cylinder  (c  =  2/3).  When  the  inner 
cylinder  was  rotated,  various  patterns  of  skewed  cells  resulted. 

Rao  et  aK  (1985)  investigated  natural  convective  flow  patterns  in 
horizontal  cylindrical  annuli.  They  appear  to  be  the  first  group  of 
researchers  that  fully  resolved  numerically  the  counter-rotating 
cells  for  air,  which  are  experienced  at  high  Rayleigh  numbers  in  the 
narrow  type  gaps.  An  unsteady  formulation  of  the  (2-D)  energy  and 
vorticity-stream  function  equations  was  used,  and  they  solved  the 
equations  using  an  ADI  scheme  with  central-differencing  throughout. 

For  R  =  1.175  and  a  Rayleigh  number  (Ra)  of  approximately  750,000, 
they  obtained  two  counter-rotating  cells  near  the  top  of  the  half¬ 
annulus.  They  reported  a  jump  in  the  steady-state  mean  Nusselt  number 
when  the  flow  made  the  transition  from  unicellular  to  multicellular 
flow.  Numerically,  they  could  not  resolve  any  type  of  oscillatory 
flow  for  the  larger  gap  widths,  but  they  did  experimentally  verify 
the  3-D  spiral  flow  associated  with  the  intermediate  size  gaps. 

Hessami  et  al_.  (1985)  studied  natural  convection  in  a  wide 
horizontal  annulus  with  a  radius  ratio  of  R  =  11.4.  They  obtained 
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experimental  and  numerical  data  for  air,  glycerin  and  liquid  mercury 

in  the  ranges  of  .023  £  Pr  £  10,000  and  .03  <  GrD  £  3  x  10^.  The 

u- 

influence  of  variable  versus  constant  fluid  properties  was  also 
explored  numerically. 

Experimentally,  they  observed  the  regular  kidney-shaped  pattern 
for  air,  and  a  multicellular  flow  pattern  for  the  liquid  mercury  case, 
similar  to  that  described  by  Mack  and  Bishop  (1968).  Numerically,  they 
assumed  vertical  symmetry  and  used  a  basic  central -differencing  scheme, 
except  for  the  convective  terms  which  were  discretized  by  using  a 
hybrid-differencing  technique  developed  by  Spalding  (1972).  This 
hybrid  scheme  collapsed  to  first-order  accuracy  upon  exceeding  the  mesh 
Reynolds  number  constraint.  Globally,  the  heat  transfer  computations 
for  air,  mercury  and  glycerin  did  not  change  with  variation  of  fluid 
properties.  However,  local  estimates  of  the  Nusselt  number  did 
exhibit  significant  discrepancy  for  glycerin  between  the  constant  and 
variable  fluid  property  cases. 

Ozoe  et  al_.  (1985)  performed  a  3-D  numerical  analysis  of  natural 
convection  in  a  spherical  annulus  for  Pr  =  1  and  Ra  =  500.  By 
imposing  a  sinusoidal  temperature  field  on  the  outer  cylinder,  they 
were  able  to  obtain  both  symmetrical  and  unsymmetrical  cell  formation. 
They  solved  the  governing  equations  with  an  ADI  scheme,  using 
central -difference  approximations  at  all  points. 
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2.4.  Variable  Fluid  Properties 
Since  the  mathematical  model  in  this  thesis  is  based  upon  the 
Boussinesq  approximation,  the  effects  of  variable  fluid  properties  and 
viscous  dissipation  on  laminar  natural  convective  flows  should  be 
addressed.  The  reviews  in  this  section  will  consider  some  of  these 
effects  with  regard  to  various  surface  geometries. 

Sparrow  and  Gregg  (1958)  analyzed  the  influence  to  variable  fluid 
properties  on  an  isothermal  vertical  flat  plate.  They  reported  that 
laminar  free  convection  heat  transfer  under  variable  property  conditions 
could  be  accurately  computed  by  using  constant  property  results  when 
evaluated  at  an  adequate  reference  temperature,  T  .  For  air  and 
liquid  mercury,  the  film  temperature,  Tf,  given  by 

Tf  *  (Tw  +  T J/2  (2.20) 

was  valid  for  most  applications.  However,  specific  reference 
temperature  relations  were  derived  and  are  listed  below. 

For  gases: 

T  =  T  -  .38  (T,  -  T  )  (2.21a) 

r  w  '  w  °° 

and  for  liquid  mercury: 

T  =  T  -  .3  (T  -  T  )  .  (2.21b) 

r  w  v  w  v  ' 
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Gebhart  (1962)  considered  the  effects  of  viscous  dissipation  on 
natural  convective  heat  transfer  for  vertical  surfaces  subject  to  both 
isothermal  and  uniform-flux  surface  conditions.  For  such  fluids  as 
liquid  sodium,  gases  and  water,  the  standard  practice  of  neglecting 
viscous  dissipation  for  natural  convective  type  flows  was  quite  valid, 
especially  in  the  laminar  regime.  But,  he  found  that  important  viscous 
dissipation  could  result  when  the  flow  made  the  transition  to 
turbulence.  Also,  he  discussed  the  fact  that  significant  viscous 
dissipation  might  occur  for  laminar  flows  subject  to  large 
decelerations  or  high  rotative  speeds. 

Gray  and  Giorgini  (1976)  determined  that  for  gases  and  most 
liquids,  in  geometries  such  as  vertical  plates  and  horizontal  cylinders, 
the  strict  Boussinesq  approximation  was  valid  for  Rayleigh  numbers 
(based  on  the  fluid  layer  depth)  up  to  10^,  nearly  13  decades  above 
the  transition  point  for  turbulence  (provided  ^  «  1). 

Clausing  and  Kempka  (1981)  experimentally  investigated  natural 
convective  heat  transfer  from  a  vertical  isothermal  heated  surface  to 
gaseous  nitrogen.  They  found  that  variable  properties  caused 
virtually  no  influence  in  heat  transfer  rates  in  the  laminar  regime, 
whereas  dramatic  increases  were  seen  in  the  turbulent  regime. 

Hessami  et  al .  (1984)  studied  the  effects  of  variable  fluid 
properties  on  natural  convective  heat  transfer  between  horizontal 
concentric  cylinders,  for  R  =  2.6.  They  concluded  that  for  air, 
the  constant  fluid  property  assumption  was  quite  valid  (for  Rayleigh 

4 

numbers  based  on  gap  up  to  9  x  10  ),  and  could  probably  be  extended  to 


all  gases.  In  contrast,  for  glycerin,  significant  differences  in  the 
temperature  field  resulted  between  the  constant  and  variable  fluid 
property  assumptions,  although  average  values  of  the  Nusselt  number 
were  very  similar. 

Mahony  ejt  aj_.  (1985)  investigated  variable  property  effects  on 
the  laminar  natural  convection  of  air  between  horizontal  cylindrical 
annuli.  They  numerically  computed  velocity  and  temperature  profiles 
for  R  =  1.5,  2.28,  2.6  and  5.0,  with  Rayleigh  numbers  based  on  gap 

5 

width  up  to  1.8  x  10  .  They  claimed  that  the  Boussinesq  approximation 

was  valid  for  a  temperature  difference  ratio,  9  ,  of  less  than  .2, 

where  9  was  given  as 
o  3 


But  for  all  the  numerical  studies  discussed  in  the  previous  section, 

0q  was  smaller  than  .1;  hence  any  related  heat  transfer  results  were 
not  affected  by  the  constant  property  assumption.  They  also  mentioned 
that  since  relatively  low  velocities  are  encountered  in  laminar 
natural  convection,  the  variable  property  assumption  should  usually 
prove  of  little  influence  on  calculated  heat  transfer  rates. 

2.5.  Flow  Bifurcation 

Due  to  the  tendency  of  multicellular  flow  to  occur  between  narrow 
horizontal  concentric  cylinders  at  high  Rayleigh  numbers,  flow 
bifurcation  with  related  hysteresis  behavior  is  definitely  possible. 


This  section  will  briefly  touch  upon  some  articles  that  aid  in  the 
description  of  nonuniqueness  and  the  bifurcation  phenomena  associated 
with  various  flows. 

Coles  (1965)  experimentally  studied  flow  transition  between 
vertical  concentric  rotating  cylinders.  The  so-called  Taylor 
instability  resulted  when  the  inner  cylinder  achieved  some  critical 
speed.  At  certain  speeds.  Coles  was  able  to  observe  both  singly  and 
doubly-periodic  type  motion.  In  this  paper,  Coles  claims  that  "the 
property  of  nonuniqueness  is  most  vividly  demonstrated  by  the 
existence  of  a  number  of  hysteresis  loops,  in  which  the  flow  changes 
from  one  state  to  another  and  back  again  as  the  speed  of  the  inner 
cylinder  is  slowly  increased  and  then  decreased."  He  mentions  that  this 
same  kind  of  behavior  is  possible  with  the  cellular  convection  patterns 
that  occur  between  horizontal  heated  surfaces. 

In  Benjamin  (1978),  some  of  the  basic  theory  associated  with 
bifurcation  phenomena  in  steady  flows  is  described.  First,  he  points 
out  that  an  instability  phenomenon  probably  exists  if  a  precise 
critical  value  of  some  parameter  can  be  related  to  the  onset  of 
cellular  motion.  Then,  a  secondary  mode  of  motion  is  usually  realized 
after  the  primary  flow  becomes  unstable  due  to  some  type  of  disruptive 
instability.  He  also  stated  that  the  Rayleigh  number  (in  the  Benard 
problem)  played  a  role  similar  to  that  corresponding  to  the  Reynolds 
number  in  the  Taylor  problem. 

Benjamin  and  Hull  in  (1982)  observed  fifteen  different  kinds  of 
steady  multicellular  flow  produced  in  a  Taylor  apparatus  with  the 


outer  wall  stationary.  They  discussed  various  forms  of  bifurcation 
that  might  result  due  to  sudden  changes  in  the  flow  field  induced  by 
variations  in  the  Reynolds  number.  And,  they  noted  that  depending  on 
the  particular  cellular  mode  initiated  at  the  outset,  many  different 
paths  to  turbulence  could  be  followed. 

Nandakumar  and  Masliyah  (1982)  investigated  the  occurrence  of  dual 
solutions  in  curved  ducts  through  a  numerical  solution  of  the  Navier- 
Stokes  equations  in  a  bipolar-toroidal  coordinate  system.  The  Dean 
number  was  the  critical  parameter  used  in  this  study,  and  was  given  by: 

On  =  Re/Rc1/2  (2.23) 

where  Re  represented  the  Reynolds  number  and  Rc  the  dimensionless  radius 
of  curvature  of  the  duct.  In  addition  to  the  Dean  number,  the  shape  of 
the  duct  was  also  varied  systematically  in  order  to  study  the 
bifurcation  of  a  two-vortex  solution  into  a  two  and  four-vortex 
solution.  They  found  that  flow  bifurcation  was  possible  irrespective 
of  the  shape  of  the  tube,  but  it  was  much  easier  to  obtain  a  dual 
solution  when  the  outer  surface  of  the  duct  was  nearly  flat. 

Cliffe  (1983)  numerically  studied  the  flow  in  a  Taylor  apparatus 
where  the  length  of  the  annulus  was  shortened  so  that  only  one  or 
two  Taylor  cells  would  result.  He  solved  the  Navier-Stokes  equations 
with  a  finite-element  method  and  then,  applied  the  methods  of 
bifurcation  theory  (see  Keller,  1977)  to  obtain  multiple  solutions  of 
the  equations  as  the  Reynolds  number  and  aspect  ratio  varied.  At  a 
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Reynolds  number  of  175,  he  obtained  three  distinct  flows:  a  stable 
two-cell  flow,  an  unstable  asymmetric  flow,  and  a  stable  single-cell 
flow.  His  numerical  method  was  able  to  resolve  the  extremely  delicate 
hysteresis  effect,  and  he  claimed  that  his  method  was  powerful  enough 
to  capture  the  more  complicated  flows  observed  in  the  Taylor  experiment 
for  moderate  aspect  ratios. 

Nandakumar  et  a]_.  (1985)  studied  laminar  mixed-convective  flow  in 
horizontal  ducts  of  rectangular,  circular  and  semicircular  cross-sections. 
In  all  cases,  dual  solutions  of  two  and  four- vortex  patterns  were 
observed.  The  governing  equations,  subject  to  the  Boussinesq 
approximation  and  an  axially  uniform  heat-flux  condition,  were  solved 
numerically  with  central -differencing  used  for  both  the  diffusive  and 
convective  terms.  For  the  case  of  the  rectangular  duct  with  a  21  x  21 
uniform  grid,  flow  hysteresis  with  respect  to  both  average  Nusselt 
number  and  friction-factor  occurred  when  the  flow  made  the  transition 
from  a  two  to  a  four-vortex  steady  solution.  The  fluid  in  this  case 
was  air,  Pr  =  .73. 

Kolodner  et  aj_.  (1986)  experimentally  studied  the  flow  patterns 
associated  with  Rayleigh-Benard  convection  in  rectangular  containers 
having  an  intermediate  aspect  ratio  of  about  10  to  5,  for  Prandtl 
numbers  between  2  and  20.  In  their  experiments,  they  observed  2-D 
skewed-varicose  and  knot  type  instabilities,  which  were  found  to 
trigger  successive  transitions  between  time-independent  flow  patterns. 

In  the  larger  Rayleigh  number  regime,  for  Pr  <  10,  the  flow 
instabilities  appeared  to  have  an  intrinsic  oscillatory-like  time 
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dependence. 

Also  related  to  flow  bifurcation  is  the  phenomena  of  strange 
attractors.  In  some  systems  that  experience  a  bifurcation  by  way  of 
a  flow  instability,  an  ordered  route  to  chaos  has  been  reported  (see 
Ruelle  and  Takens,  1971;  Newhouse  et  al_. ,  1978;  Giglio  et  al . ,  1981; 
Brandstater  et  aj_. ,  1983;  Grebogi  et_  al_. ,  1983;  and  Guckenheimer , 
1986).  A  typical  sequence  of  events  is  as  follows.  First,  the  system 
behaves  in  a  time-periodic  manner  after  the  onset  of  the  initial 
instability.  Then,  upon  further  increase  of  a  system  parameter  (such 
as  Reynolds  or  Rayleigh  number),  a  cycle  of  periodic-doubling  is 
usually  observed,  until  finally,  chaotic  behavior  sets  in  via  small- 
scale  spectral -broadening. 

Some  of  these  strange  but  ordered  characteristics  seem  to  have 
been  numerically  predicted  in  this  present  study.  A  description  of 
the  numerical  technique  and  related  results  will  be  given  in 
Chapters  5  and  6,  respectively. 

2.6.  Natural  Convection  in  Vertical  Slots 

Based  upon  the  work  set  forth  in  this  thesis,  it  appears  that  a 
multicellular  instability  may  occur  in  the  vertical  sections  of  very 
narrow  horizontal  cylindrical  annuli.  Because  of  this  likelihood, 
several  papers  dealing  with  laminar  natural  convection  in  vertical 
slots  will  be  reviewed  in  order  to  shed  more  light  on  the 
manifestation  of  this  unique  type  of  hydrodynamic  instability. 
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Elder  (1965)  performed  an  experimental  study  of  natural  convection 
for  a  liquid  (Pr  =  1,000)  in  a  vertical  slot  with  isothermal  walls 
(the  left  wall  being  hottest).  Unicellular  flows  resulted  when  the 

5 

Rayleigh  number  based  on  gap  width  was  less  than  10  .  But  for  the 
larger  aspect  ratio  (narrow  gap)  slots,  with  Rayleigh  numbers  slightly 
greater  than  10  ,  a  steady  'cats-eye'  patterned  secondary  flow  became 
superimposed  on  the  basic  unicellular  flow.  And,  upon  further 
increasing  Rayleigh  number  to  above  10^,  Elder  (1965)  reported  the 
emergence  of  a  tertiary  flow  with  counter-rotating  type  cells. 

Elder  (1966)  numerically  solved  the  same  problem  as  described 
above,  except  in  this  case,  he  considered  only  the  moderate-size 
vertical  slots.  He  was  able  to  duplicate  the  basic  flow  field 
obtained  in  his  experiments,  but  was  not  able  to  resolve  the  secondary 
flow  that  was  present  in  his  previous  work.  Elder  mentioned  that  the 
nonlinear  terms  began  to  dominate  the  motion  as  Rayleigh  number 
increased,  and  he  was  able  to  numerically  show  the  development  of  the 
boundary-layers  and  the  fully-developed  boundary-layer  flow. 

Vest  and  Arpaci  (1969)  analytically  investigated  the  stability 
of  natural  convection  in  a  narrow  vertical  slot.  By  using  linearized 
hydrodynamic  stability  theory,  they  were  able  to  obtain  a  neutral 
stability  curve  for  the  conduction-dominated  flow  regime.  For  .01  <  Pr 
<_  10,  the  variation  of  the  critical  Grashof  number  was  found  to  be  less 
than  .7  percent.  Hence,  for  Prandtl  numbers  in  this  range,  they 
determined  a  single  critical  Grashof  number  of  7,880  at  a  wavenumber 
of  2.65.  Also,  their  analytically  obtained  stream  function  plots  for 
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the  secondary  flow  qualitatively  agreed  with  related  flow  pictures 
obtained  experimentally. 

Thomas  and  Vahl  Davis  (1970)  numerically  studied  natural  convection 
between  vertical  cylindrical  annuli,  with  the  inner  cylinder  being  at 
a  higher  temperature  than  the  outer.  (Note  that  a  similar  study  with 
more  basic  results  can  be  found  in  Vahl  Davis  and  Thomas,  1969.) 

They  solved  the  unsteady  vorticity  and  energy  equations  using  an  ADI 
scheme.  For  an  aspect  ratio  of  H  =  25  (length  of  annulus  to  gap 
width),  and  a  Rayleigh  number  of  22,500  (based  on  gap  width),  they 
observed  an  unsteady  type  of  multicellular  flow  for  Pr  =  1.0.  This 
phenomena  was  similar  to  that  reported  by  Elder  (1965),  except  that 
Elder's  secondary  flow  maintain'd  the  steady-state  condition. 

Korpela  et  al_.  (1973)  used  linear  stability  theory  to  examine 
the  stability  of  the  conduction  regime  for  natural  convection  in  a 
vertical  slot.  For  Pr  <  12.7,  they  claimed  that  the  instability  set 
in  as  horizontal  cells.  They  observed  a  critical  Grashof  number  of 
7,932,  at  a  wavenumber  of  2.65,  for  Pr  =  0.  This  particular  type  of 
instability  was  thought  to  be  hydrodynamic  in  origin,  resulting  from 
the  vorticity  distribution  of  the  base  flow. 

Korpela  (1974)  studied  a  problem  similar  to  that  described  above, 
except  in  this  case,  he  assumed  that  the  narrow  slot  was  maintained 
at  an  angle  6  with  the  vertical.  For  Pr  <  12.7,  he  found  that  the 
instability  set  in  as  transverse  travelling  waves  for  small  angles  of 
inclination,  and  that  longitudinal  cells  formed  as  6  reached  a  certain 
value.  In  the  range  of  .24  <_  Pr  <  12.7,  he  determined  that  the 
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instability  would  lead  to  horizontal  cells  for  angles  close  to  the 
vertical,  whereas  longitudinal  cells  would  result  as  the  slot  was 
further  inclined.  For  Pr  <  .24,  he  claimed  that  only  horizontal 
cells  were  possible  and  that  the  stability  of  the  flow  was  mainly  a 
function  of  Pr  tan  5. 

Pepper  and  Harris  (1977)  numerically  obtained  2-D  natural 

convective  flow  patterns  in  rectangular  and  annular  vertical  cavities. 

The  energy  and  vorticity  equations  were  written  in  divergence  form  and 

were  solved  using  central -differencing  with  a  strongly  implicit 

proo.-f  re.  ror  the  rectangular  slot  with  an  aspect  ratio  of  10,  for 

Pr  =  1,000,  they  obtained  a  weak  multicellular  flow  pattern  at  a 
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Rayleigh  number  (based  on  gap  width)  of  approximately  5  x  10  . 

Seki  et  aj_.  (1978)  performed  an  in-depth  experimental  analysis  of 
natural  convection  in  narrow  vertical  rectangular  cavities.  They 
considered  transformer  oil,  water  and  glycerin  as  the  working  fluids. 
All  three  fluids  yielded  a  multicellular  type  of  secondary  flow  at  a 
certain  Rayleigh  number,  Ra^,  based  on  the  height  of  the  slot.  For 
oil,  as  the  temperature  difference  was  increased,  they  noticed  a 
tertiary  flow  with  counter-rotating  cells,  until  finally,  the  flow 
near  the  upper  region  of  the  hot  wall  became  unsteady  and  turbulent. 

For  the  oil,  for  an  aspect  ratio  of  15,  the  secondary  motion 

g 

began  at  Ra^  =  6  x  10  and  the  transition  to  turbulence  took  place  at 
about  Ra^  =  1.5  x  10^.  They  concluded  that  the  flow  field  more  easily 
shifted  from  laminar  to  transitional  flow  as  the  Prandtl  number 
decreased  and  the  cavity  width  increased.  Another  interesting 
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description  of  these  secondary  and  tertiary  motions  can  be  found  in 
Vahl  Davis  and  Mallinson  (1975). 

Using  linear  theory,  Choi  and  Korpela  (1980)  studied  the  stability 
of  the  conduction  regime  for  natural  convection  in  a  vertical  annulus. 
They  found  that  for  all  Prandtl  numbers,  the  instability  set  in  as 
an  upward  travelling  wave.  Hence,  stationary  cells  were  no  longer 
possible  as  with  the  vertical  slot  geometry.  For  low  Prandtl  numbers, 
the  larger  the  curvature,  the  more  stable  the  flow,  while  the  reverse 
was  true  for  the  higher  Prandtl  number  fluids.  From  experimental  data, 
the  measured  wavelength  of  the  cells  was  in  good  agreement  with  their 
linear  analysis. 

Orszag  and  Kells  (1980)  studied  the  role  of  two-  and  three- 
dimensional  finite-amplitude  disturbances  in  the  breakdown  of  plane 
Poiseuille  and  plane  Couette  flows.  To  determine  the  evolution  of 
these  disturbances,  they  solved  the  3-D  time-dependent  Navier-Stokes 
equations  using  spectral  methods  with  Fourier  and  Chebyshev  polynomial 
series.  They  claimed  that  the  3-D  finite-amplitude  effects  produced 
strong  inflexional  velocity  profiles  that  eventually  caused  the 
transition  to  turbulence;  whereas  the  2-D  disturbances  proved  to  be 
much  less  destabilizing,  and  seemed  powerless  for  Reynolds  numbers 
below  3,000.  The  contour  velocity  plots  due  to  the  3-D  disturbances 
evolved  into  a  'cats-eye'  pattern  very  similar  to  that  observed  by 
Lee  and  Korpela  (1983)  and  Elder  (1965)  in  their  study  of  multicellular 
natural  convective  flow  in  a  vertical  slot. 


Lee  and  Korpela  (1983)  succeeded  in  numerically  resolving 
multicellular  natural  convective  flow  in  narrow  vertical  slots.  They 
attributed  their  success  to  the  4th-order  differencing  method  of 
Arakawa  (1966),  which  was  used  in  approximating  the  nonlinear 
convective  terms.  The  important  buoyancy  term  for  the  vertical  slot 
geometry  was  represented  by  the  single  component,  3T/3x,  in  the 
vorticity-transport  equation.  This  resembled  the  buoyancy  term 
obtained  in  the  boundary-layer  equations  of  Chapter  4  of  this  thesis, 
where  the  term  reduced  to  simply  3T/3r  at  ^  =  90°  (see  Figure  3.1  and 
Eq.  (4.25b)).  They  were  able  to  obtain  steady  multicellular  flow  for 
Prandtl  numbers  ranging  from  zero  to  1,000.  However,  as  the  Prandtl 
number  increased,  the  aspect  ratio  (H)  had  to  be  increased 
significantly  in  order  to  trigger  the  multicells.  For  Pr  =  0,  and 
H  =  15,  they  obtained  a  secondary-flow  transition  to  six  cells  at  a 
Grashof  number  (based  on  gap  width)  of  8,000.  This  number  agreed 
quite  well  with  their  analytical  prediction  of  7,932  derived  in 
Korpela  et  al_.  (1973).  The  six  cells  formed  in  the  vertical  slot 
appeared  fairly  constant  in  strength.  They  were  not  sure  whether  a 
further  transition  to  periodic  flow  was  possible  for  this  2-D  flow  or 
whether  the  next  physically  important  structure  was  a  3-D  steady  or 
periodic-type  flow. 


2.7.  Concluding  Remarks 

This  literature  review  discusses  the  various  assumptions  and 
solution  methods  employed  by  researchers  in  studying  natural  convection 


between  horizontal  concentric  cylinders.  Some  contradictory 
conclusions  and  results  have  been  highlighted.  Of  particular 
importance,  experimenters  have  verified  that  for  air  at  high  Rayleigh 
numbers,  a  multicellular  flow  is  possible  near  the  top  of  narrow 
horizontal  cylindrical  annuli.  However,  not  all  numerical  investigators 
confirmed  this  fact.  It  seemed  that  the  transition  to  multi  cells  could 
not  be  captured  due  to  either  first-order  upwind  differencing  of  the 
convective  terms,  or  a  lack  of  numerical  stability  experienced  when 
approaching  high  Rayleigh  numbers.  Thus,  it  appears  that  at  least 
second-order  accuracy  is  a  prerequisite  for  resolving  this  multicellular 
type  of  instability.  In  addition,  for  this  geometry  there  appears 
to  be  a  void  in  the  research  regarding  the  possibility  of  related 
hysteresis  behavior  associated  with  the  transition  to  multicells. 

Also  from  this  literature  review,  it  can  be  seen  that  the 
potential  multicellular  flow  in  the  vertical  portions  of  narrow 
horizontal  annuli,  as  described  in  this  research  effort,  has  not  been 
examined. 


3.  MATHEMATICAL  ANALYSIS 


In  this  study,  the  Boussinesq  approximated  Navier-Stokes  equations 
and  the  viscous-dissipation  neglected  thermal -energy  equation  will  be 
used  to  determine  the  buoyancy-induced  steady  or  unsteady  flow  fields 
between  horizontal  isothermal  concentric  cylinders.  The  Boussinesq 
approximation  states  that  density  perturbations  are  only  felt  in  the 
body  force  terms  and  can  be  neglected  in  the  acceleration  and  viscous 
components.  Furthermore,  the  vorticity-stream  function  formulation  of 
the  Navier-Stokes  equations  will  be  adopted.  Using  this  approach 
completely  eliminates  the  pressure  gradient  terms  and  automatically 
satisfies  the  conservation  of  mass  principle.  An  added  advantage  is 
that  stream  function  and  vorticity  contours  are  well -suited  for 
visualizing  and  analyzing  the  flow. 

3.1.  The  Physical  Model 

a.  The  model  is  unsteady  and  two-dimensional  (see  Figure  3.1). 

b.  Initially,  the  fluid  is  at  rest. 

c.  The  cylinders  are  assumed  horizontal  and  isothermal,  with  the 
inner  cylinder  temperature  exceeding  the  outer. 

d.  Laminar  fluid  motion  is  induced  by  buoyancy  effects.  The 
fluid  is  Newtonian. 

e.  All  material  properties  are  assumed  constant.  Density 


variations  are  allowed  to  occur  via  the  Boussinesq 
approximation. 


3.2.  The  Complete  Set  of  Governing  Equations 
Before  formally  deriving  the  dimensional  and  nondimensional  forms 
of  the  governing  equations,  the  complete  system  of  equations  will  be 
given  below.  This  is  done  to  provide  for  easy  access  when  the  equations 
are  numerically  solved  or  referred  to  in  subsequent  chapters  of  this 
thesis.  The  dimensionless  system  of  coupled  partial  differential 
equations  is 

Thermal -energy: 


r2  II  4.  4.  1\-1  /3£  U  9f  £]\  -  JL  rl!l 

b  3t  1  G'  l3r  "  dip  dr'  '  Pr  ^2 


+  (r  +  If  +  (r  +  s}’2  f?} 


(3.1a) 


Vorticity : 


p_  3w  /  lrl  / 3f  3w  3f  3wn 
Pr  {G  (r  +  q)  < aF  3^7  "  3^7 


=  Pr  r&L  +  ,  lrl  aw  ,  K-2  3^w  } 

rr  \  2  lr  GJ  3r  +  ir  G  .,2  1 

3r  dip 

+  G(Ra)  {sir  *  |I  + 31  ) 


(r  +  i) 


(3.1b) 


Stream  function: 


4  Mr  +  I)’’  |f  Mr  +  I)’2  0  ■  G2* 

3r  dip 


(3.1c) 


while  the  related  boundary  conditions  are 


at  r  =  0  (inner  cylinder): 


T  =  1 

2 

u,  -  1  9  f 

w  _  — 2 
3r 


f  =  0 


(3. Id) 


at  r  =  1  (outer  cylinder) 


T  =  0 


G  3r 


f  =  0  . 


(3. le) 


The  derivation  and  explanation  of  important  terms  and  symbols,  together 
with  the  boundary  and  initial  conditions,  follow  in  the  next  two 
sections. 


3.3.  The  Dimensional  Formulation 
3.3.1.  Governing  equations  in  primitive  variables 


Four  governing  equations  are  needed  in  order  to  describe  the 
transport  of  mass,  momentum  and  energy  between  the  horizontal 
concentric  cylinder  geometry.  These  equations  are  written  in  terms  of 
primitive  variables  using  polar  coordinates  ( r ,^ )  (see  Figure  3.1). 

The  bar  over  the  variable  signifies  a  dimensional  quantity.  The  radial 
and  angular  coordinates  are  given  by  r  and  iJj,  respectively.  T  is  the 


temperature  and  t  denotes  time.  The  radial  velocity  is  represented  by 
u  and  the  angular  or  tangential  velocity  is  given  by  v.  P  is  the 
pressure  term,  while  v  and  a  are  the  momentum  and  thermal  diffusi vities , 
respectively.  Lastly,  p  is  the  density  and  F  is  the  body  force  term 
in  the  equations  of  motion.  The  2-D  governing  equations  in  primitive 
variables  are 

Equation  of  continuity: 

iy.  +  y.  +  lil=0  (3.2a) 

3r  r  r 

Equations  of  motion: 

2.  d.3£  +  v(v2-  i  2  2i,  +  f  (3.2b) 

Dt  r  0  3r  r  r  3'U 

4  +  S*.  4  3£  t  v(v2;  .  J  t  2  li)  +  f  (3.2c) 

Ot  r  pr  3ip  r  r  dip  * 

Equation  of  thermal -energy: 

=  a  V2  T  (3. 2d) 

Dt 

where  D/Dt,  the  substantial  or  material  derivative,  represents  a 
combination  of  the  local  and  convective  changes  of  a  particular  property 
(usually  associated  with  the  acceleration  terms),  and  is  given  by 


Also,  v  (the  2-D  cylindrical  Laplacian)  is  given  by 


Z  _  32  ,  1  9  .1  3 

+  ~  ~  +  ^2  ~?i  • 

3r  r  3r  r  3^ 


(3.3) 


“2  * 

The  centrifugal  acceleration  term,  v  /r,  in  Eq.  (3.2b),  represents 
the  effective  force  per  unit  mass  and  volume  in  the  radial  direction, 
resulting  from  fluid  motion  in  the  angular  direction.  Whereas,  uv/r  in 
Eq.  (3.2c)  designates  the  Coriolis  acceleration,  and  represents  the 
effective  force  (per  unit  mass  and  volume)  in  the  angular  direction 
when  fluid  motion  occurs  in  both  the  radial  and  angular  directions. 

It  should  be  noted  that  these  two  terms  arise  automatically  upon 
transformation  from  rectangular  to  cylindrical  coordinates  (Bird  et  al . 
1960).  The  F  terms  in  Eqs.  (3.2b)  and  (3.2c)  are  due  to  the  buoyancy 
force.  They  can  be  evaluated  with  the  aid  of  the  Boussinesq 
approximation,  which  allows  density  to  be  variable  only  in  the  body 
force  terms  of  the  momentum  equation.  The  corresponding  change  in 
fluid  density  with  temperature  is  related  by 


p'  =  p  -  po  =  -pQe  (T  -  tq) 

where  pQ,  Tq  represent  property  values  evaluated  at  some  known 
reference  state  (e.g.,  the  initial  density  and  temperature  of  the 
fluid  within  the  annulus),  and  8  is  the  coefficient  of  thermal 
expansivity,  which  becomes  the  negative  reciprocal  of  absolute 
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temperature  for  an  ideal  gas.  One  can  regard  the  expression  for  p  as 

the  first  two  terms  of  a  Taylor-series  expansion  about  T  .  The  net 

i  body  force  term,  g  j,  can  now  be  broken  down  into  its  radial  and 

!  po 

angular  components: 

i 

i 

-gS(T  -  TQ) j  =  -ge(T  -  T0)cosi|)  er 

+  9B(T  -  TQ)simj;  e^  (3.4) 

where  g  is  the  acceleration  of  gravity. 

The  Boussinesq  approximation  is  normally  considered  quite  accurate 
when  dealing  with  laminar-type  flows,  especially  for  Prandtl  numbers 
ranging  from  near  zero  to  one  and  for  small  AT  (AT  ~  0  (15°C)  or  less) 
(Gray  and  Giorgini,  1976;  Hessami  et  al_. ,  1984;  Mahony  et  ,  1985). 

As  a  final  note,  the  viscous  dissipation  term  has  not  been  included 
in  the  thermal -energy  equation.  This  term  is  usually  only  significant 
in  high-speed  flow  applications,  or  for  laminar  flows  subject  to  large 
decelerations  (Gebhart,  1962). 


3.3.2.  Governing  equations  using  the  stream  function-vorticity  approach 
The  stream  function  (f)  is  defined  in  such  a  way  so  as  to 
automatical ly  satisfy  the  continuity  equation.  It  is  valid  for  all 
two-dimensional  flows,  both  rotational  and  irrotational ,  and  can  be 
defined  as: 


u  -  4  ^  and 


v  = 


3f 


( 


3.5 


The  number  of  dependent  velocity  variables  has  now  been  reduced  by  one. 
Physically,  for  2-D  incompressible  flows,  the  volumetric  flow  rate 
between  two  streamlines  is  given  by  the  difference  in  f. 

To  further  simplify  the  remaining  governing  equations,  the 
vorticity-transport  equation  will  be  derived  by  taking  the  curl  of  the 
Navier-Stokes  equations.  To  begin,  the  equations  of  motion  for  a 
Newtonian  fluid  of  constant  density  and  viscosity  will  be  written  in 
vector  form  (which  is  valid  for  any  coordinate  system): 

~r  +  (v  •  V)v  =  +  vV2  v  +  F0  .  (3.6) 

o  t  P  D 

Note  that  V  =  +  -  —  e,,  in  2-D  cylindrical  coordinates,  and  v 

3r  r  r  3if>  * 

represents  the  velocity  vector.  Next,  the  convective  term  in  Eq.  (3.6) 
is  expanded  using  the  following  vector  identity: 

(v  •  V)v  =  1  V(v  •  v)  -  v  x  (V  x  v)  . 

Then,  since  vorticity  (which  physically  represents  the  rotation  of 
an  infinitesimal  fluid  particle)  is  given  by 

w  =  V  x  v  , 

and  by  utilizing  other  appropriate  vector  identities  along  with  the 
continuity  equation,  one  obtains  the  following  result: 
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!y  +  (v  •  V)w  =  (w  •  7)v  +  W2w  +  V  x  Fb  .  (3.7) 

The  left  side  on  Eq.  (3.7)  describes  the  total  rate  change  of  particle 

vorticity.  (w  •  V)v,  which  represents  the  rate  of  stretching  of 

vortex  lines,  is  identically  zero  for  2-D  flows  since  the  vorticity 

vector  in  the  2-D  case  is  always  perpendicular  to  the  plane  of  flow. 

2- 

vV  w  signifies  the  net  rate  of  vorticity  diffusion  due  to  viscous 
effects.  The  last  term,  V  x  Fg,  represents  the  rate  of  internal 
vorticity  generation  due  to  body  forces,  which  result  from  significant 
density  perturbations  in  natural  convective  flows  (Panton,  1984). 

Note  the  fact  that  the  pressure  term  does  not  appear  explicitly  in  the 
vorticity  equation.  Thus,  the  vorticity  and  stream  function  (velocity) 
fields  can  be  determined  without  any  prior  knowledge  of  the  pressure 
distribution  (Currie,  1974). 

Using  the  result  of  Eq.  (3.4)  to  evaluate  V  x  Fg,  and  upon 
substituting  the  stream  function  for  the  velocity  terms,  the  final  form 
of  the  vortici ty-transport  equation  is  obtained,  namely: 

%  +  l-  =  gg(sint|j  4  +  4)  +  W2w  (3.8a) 

3t  r  3(r,40  3r  r  3^ 


and  likewise,  the  thermal -energy  and  stream  function  equations 
simplify  to: 


(3.8b) 


where  V  is  the  same  as  in  Eq.  (3.3),  and  the  Jacobian,  3(P,Q)/3(x,y) , 


is  defined  as 


3(P,Q)  =  3P  3Q.  _  3P  li 

3(x,y)  3x  3y  3y  3x  ‘ 


Equations  (3.8a,  b,  and  c)  are  the  three  coupled  governing  equations 
that  describe  the  dependent  variables  w,  T,  and  f  for  the  horizontal 
annulus. 


3.3.3.  Boundary  conditions 

The  energy,  vorticity,  and  stream  function  equations  are  all 
second-order  partial  differential  equations.  Therefore,  two  boundary 
conditions  (in  each  spatial  coordinate)  for  each  equation  are  needed 
to  properly  define  a  well-posed  mathematical  problem.  For 
temperature,  isothermal  cylinders  are  assumed,  with  the  inner 
cylinder  being  hotter  than  the  outer  one.  Thus, 


T  =  T^  at  r  =  a 


T  =  Tq  at  r  =  b, 


(3.9a) 


where  T^  is  greater  than  T  . 

Since  both  u  and  v  are  identically  zero  at  the  stationary  cylinder 
walls  due  to  the  no-slip  condition,  the  vorticity  boundary  condition 


J'fcJU'LM 


becomes : 


at  r  =  a  and  r  =  b 


(3.9b) 


Again,  because  of  no-slip  at  the  walls,  the  stream  function  must  remain 


constant  and  can  be  taken  as: 


f  =  0  at  r  =  a  and  r  =  b  . 


(3.9c) 


Since  the  complete  annular  cylindrical  space  (0-2tt)  is  being  analyzed, 
the  following  boundary  conditions  (in  ^)  are  needed  because  of 
continuity  at  zero  and  2ir  radians: 


^rr 


1£  =  3£ 

o  d'P  2tt 


for  any  given  r,  a  <_  r  <  b  ,  (3.9d) 


where  <p  takes  on  values  of  either  temperature,  vorticity,  or  stream 
function. 

For  certain  cases  (as  will  be  described  in  the  numerical  analysis 
of  Chapter  5),  syrrenetry  about  the  vertical  axis  will  be  assumed.  In 
these  cases,  the  condition  defined  by  Eq.  (3.9d)  will  be  replaced  by 
the  following  symmetrical  boundary  conditions: 


f  =  w  =  0  and  —  =  0  at  =  0  and  tt  radians 


(3.10) 
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3.3.4.  Initial  conditions 

Since  initially  there  is  no  fluid  motion,  the  initial  conditions 
for  this  problem  are  straightforward.  Thus, 

f  =  w  =  0  and  T  =  Tq  (3.11a) 

throughout  the  annulus,  except  at  the  walls  where: 

f  =  T.  and  T  =  Tq  (3.11b) 

at  r  =  a  and  r  =  b,  respectively. 

Prescribed  initial  conditions,  other  than  "zero"  fluid  motion, 
will  be  used  in  this  study  when  checking  for  hysteresis  loops.  These 
conditions  will  be  briefly  discussed  after  the  nondimensional  formulation 
is  derived  in  the  next  section. 

3.4.  The  Dimensionless  Formulation 

The  advantages  of  casting  equations  in  dimensionless  format  are 
that  first,  it  helps  transform  the  mathematical  or  numerical  results 
into  a  simpler  (normalized)  form,  thus  allowing  for  better  graphical 
interpretation.  Second,  measurement  scales  are  no  longer  an  intrinsic 
part  of  the  physical  quantities,  so  that  any  laws  governing  physical 
variables  are  valid  tor  all  diMerent  measurement  systems  (Panton, 

1984,.  Most  important1,.  «rvr  a  ;-r-->ti>r  ’  $  nondimensional  i  zed,  fewer 
variables  are  js*>d  .vv,:  •’  >  .  ■  :>  •  :  ■  n i ess  jroups  needed  to 

specify  a  par*;  ./ v  ,  * .  .  T  ’  (zr*  h. 


3.4.1.  Coordinate  transformation 


The  radial  coordinate  is  nondimensional ized  so  that  the  outer 
boundary  at  r  =  b  is  transformed  into  the  unit  circle  r  =  1  (see 
Figure  3.2).  The  inner  boundary,  r  =  a,  is  transformed  into  the  pole, 
r  =  0.  The  following  is  then  obtained: 


r  =  £-—•  ,  <P  =  4>  (coordinate  (3.12a)  ; 

transformation)  j 

i 

i 

and  the  other  independent  variable,  time,  is  scaled  as  ! 


t  =  t  ^  .  (3.12b) 


{ 

3.4.2.  Governing  equations  j 

t 

The  remaining  dependent  variables  are  now  scaled  to  produce  the  1 

i 

following  dimensionless  set  of  dependent  variables:  - 


(temperature)  (3.12c)  ! 

t 


(vorticity)  (3.7  2d )  J 

< 


(stream  function)  .  (3.12e)  j 

| 

i 

Equations  (3.12)  can  now  be  substituted  into  Eqs.  (3.8a-c),  j 

1 

yielding  the  following  nondimensional  governing  equations:  j 

i 

| 

i 

\ 

\ 

) 
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Pr  {g2  ||  ♦  (r  ♦  If'  0^}  >  ■  f  •  v|w 


+  G  •  Ra  {sin^  + 


9T  .  cos^  9T 


(r 


lx  W 


r2  9T  .  /  ,  1\_1  9(f>T)  _  1  /n2 

G  at  +  (r  G}  w:$J~  ^  (  2  T) 


v2  f  =  G2w 


(3.13a) 


(3.13b) 


(3.13c) 


where 


(r 


+  l)'1  1_ 

G  9r 


tr  ♦  I)’2 


9_ 

9^ 


2 


and  again. 


9(P,Q)  _  9P  3Q  9P  92. 
9(x,y)  9x  3y  "  9y  3x 


Also, 


(Prandtl  number) 


Ra 


9Ba3(Ti  -  T0) 

v  a 


(Rayleigh  number 
based  on  the  inner 
radius,  a) 


(3.13d) 

(3.13e) 


and  G  =  — "  a  (gap  number)  .  (3.1 3f ) 

G  was  chosen  such  that  as  the  inner  radius  (a)  approached  the  outer 
radius  (b),  the  gap  number  would  tend  to  zero. 

\ 

i 

] 


Therefore,  the  three  dimensionless  groups  brought  forth  from  this 


analysis  are  G,  Pr,  and  Ra.  These  three  variables  will  be  used  to 
simulate  the  various  flow  conditions  and  geometries  under  consideration 
The  Rayleigh  number  alone  is  normally  considered  the  critical  parameter 
which  predicts  the  onset  of  thermal  and/or  hydrodynamic  instability 
(variations  of  this  will  be  discussed  in  Chapter  4),  and  can  be 
interpreted  as  the  ratio  of  the  destabilizing  buoyant  forces  to  the 
stabilizing  viscous  forces.  The  Prandtl  number  characterizes  the 
fluid,  and  the  gap  number  identifies  the  geometry. 

In  the  formulation  of  Eqs.  (3.13),  the  variables  w,  f,  and  t  were 
nondimensionalized  with  v  instead  of  a.  This  was  done  in  the  hope  that 
for  Pr  <  1 ,  the  numerical  stability  of  the  problem  would  be  enhanced. 
For  these  low  range  Prandtl  numbers,  investigated  in  this  study,  the 
flow  field  should  be  thermal -diffusion  dominated.  This  is  consistent 
with  the  thermal -energy  equation  (Eq.  3.13b);  that  is,  as  Pr  decreases 
from  one,  the  diffusion-side  of  Eq.  (3.13b)  increases.  Representing 
the  equations  in  this  nondimensional  form  better  supports  diagonal 
dominance  when  the  equations  are  solved  numerically,  since  less 
emphasis  is  given  to  the  left-hand  side  (off-diagonal  contributing) 
nonlinear  terms. 


The  Jacobian  terms  in  Eqs.  (3.13a  and  b)  can  also  be  written  in 
terms  of  the  velocity  variables,  where  u  and  v  are: 


— —  ■)  =  u  (-)  (radial  velocity) 


(3.14a 


(tangential  velocity). 


With  these  substitutions,  the  dimensionless  governing  equations  in 
(3.13)  become: 


Energy: 


r2  9T  .  „r  9T  .  vG  9T 
G  «  +  uG^+(r  + 1,» 


,  1_  r  afl  .  1  £L  .  1  32t 

Pr  ar2  (r  +  J-)  ar  (r  +  ^)2  9^2 


Vorticity: 


Pr  {G2  |w  t  uG  |S  +  |W  , 

-  pr  f3  w  ,  1  9w  1  3^w  i 

{a7  (r  +  l)9r  (r  7  l)2  a*2  } 


+  G(Ra)  ^in^fI  +  -^-|I} 

(r  +  3) 


Stream  function: 


92f  1  3f  .  1  52f 

9?  (r  +  i)  9r  (r  +  Jr)2  9>P2 


Equations  (3.15a-c)  were  written  in  the  above  form  so  as  to  facilitate 
the  asymptotic  analysis  that  will  follow  in  Chapter  4,  in  which  the 
double-limit  of  Ra  -*  00  and  G  0  is  examined. 


3.4.3.  Boundary  conditions 

The  nondimensional  boundary  conditions  for  this  problem  are 
obtained  by  substituting  Eqs.  (3.12)  into  Eqs.  (3.9).  They  are  given 
in  Eqs.  (3. Id)  and  (3.1e)  for  r  =  0  and  1,  respectively. 

The  continuity  boundary  condition  applies  at  zero  and  2ir  radians, 
and  is  similar  to  that  described  in  Eq.  (3.9d),  except  that 
dimensional  variables  are  now  nondimensional ized.  Likewise,  if 
vertical  symmetry  is  assumed,  one  gets: 


f  =  w  =  0,  ^  =  0  at  ^  =  0  and  it  radians  . 


(3.15 


3.4.4.  Initial  conditions 

The  dimensionless  initial  conditions  are  essentially  the  same  as 
those  listed  in  Eqs.  (3.11).  That  is,  for  the  entire  annulus. 


T  =  f  =  w  =  0  , 


(3.16 


except  at  the  isothermal  walls,  where: 


T  =  1  at  r  =  0 


T  =  0  at  r  =  1 


(3.16 


■»  A 


64 


Again,  these  initial  conditions  arise  by  assuming  that  initially  the 
fluid  is  motionless,  or  Ra  =  0.  Thus,  the  problem  initially  reduces 
to  the  1-D  steady-state  conduction  case  for  a  cylindrical  annulus. 
However,  when  numerically  searching  for  transitional  Rayleigh  numbers 
and  hysteresis  behavior,  initial  conditions  other  than  "zero"  are 
used.  These  procedures  will  be  outlined  in  detail  in  Chapter  5. 


4.  ASYMPTOTIC  ANALYSIS 


In  this  chapter,  a  high  Rayleigh  number/small-gap  asymptotic 
expansion  theory  is  described  in  which  the  2-D  Navier-Stokes  (N-S)^ 
equations  collapse  into  cartesian-like  boundary-layer  equations. 

The  purpose  of  this  analysis  was  to  gain  a  deeper  understanding 
of  the  transitional  tendency  toward  multicellular  flow  in  narrow  gaps, 
especially  with  regard  to  the  mechanism  that  triggers  the  multicells. 

To  understand  this  mechanism,  key  points  had  to  be  resolved:  (1)  does 
the  source  of  instability  originate  from  thermal  or  hydrodynamic 
effects,  and  (2)  once  determined,  do  these  effects  ever  coexist  or  does 
one  always  outweigh  the  other? 

Also  from  this  asymptotic  analysis,  results  were  obtained  which 
serve  as  convenient  checks  for  the  pretransitional  numerical  results 
generated  in  the  flow  bifurcation  study  of  this  research  effort. 

At  the  high  Rayleigh  number  limit,  the  annular  flow  field  can  be 
divided  into  an  inner  and  outer  boundary-layer  with  an  inviscid  core 
in  the  center.  In  general,  the  boundary-layers  are  noninteracting, 
which  disallows  an  inviscid  core  solution  to  be  arbitrarily  set. 

Also,  since  inviscid  flows  are  nonunique,  the  correct  core  solution  must 
be  chosen  in  order  to  avoid  a  singularity  at  separation  (Goldstein, 
1948).  Since  the  correct  inviscid  solution  is  usually  not  known  a 
priori,  one  needs  a  type  of  interactive  boundary-layer  theory  which 


Throughout  this  chapter,  the  2-D  N-S  equations  refer  to  those 
outl ined  in  Eq.  (3.15). 


"s' 


supports  separation.  A  good  way  of  devising  such  a  theory  for  this 
particular  geometry  is  to  look  for  the  point  where  the  boundary-layers 
merge.  Then,  the  inviscid  core  can  be  asymptotically  matched  with  the 
boundary-layer  to  verify  the  small-gap  boundary-layer  structure. 

Work  begins  with  the  construction  of  inviscid  core  expansions  for 
Ra  -*■  oo  and  G  0.  This  construction  is  necessary  to  facilitate  a  valid 
match  of  the  core  and  boundary-layer  tangential  velocities,  which,  in 
turn,  provides  for  the  correct  scaling  of  the  boundary-layer  variables. 

4.1.  The  Inviscid  Core 

For  finite  values  of  the  coordinates,  the  viscous  terms  in  the 
2-D  N-S  vorticity  equation  can  be  neglected  in  the  limit  as  Ra  -+  o°. 

This  results  in  an  inviscid  vorticity  equation  in  which  the  buoyancy 
forces  are  balanced  by  the  convective  acceleration  forces  at 
steady  state.  Since  vorticity  is  represented  as  the  curl  of  the 
velocity  field,  and  because  the  gap  scaling  has  been  taken  care  of  in 
the  nondimensionalization,  one  can  assume  that  velocity  expands  like 
vorticity,  which  gives 

w  ~  u  ~  v  ~  Ra^ 
c  c  c 

by  inspection  of  the  inviscid  vorticity  equation.  More  formally,  the 
following  expansions  can  be  substituted  into  the  inviscid  vorticity 
equation  to  provide  the  same  result  to  leading  order,  namely, 


uc 

=  <S  u 

0  0 

+ 

Vi 

+  .  •  • 

vc 

II 

o> 

O 

< 

O 

+ 

Vi 

+  •  .  . 

wc 

=  {  w 

O  0 

+ 

«,w, 

+  .  . . 

fc 

II 

.3 

O 

o 

+ 

nlfl 

+  •  .  . 

Tc 

=  Y  T 

0  0 

+ 

Y1T1 

+  •  . . 

(4.1) 
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where  SQ  ~  Ra  .  Likewise,  upon  substituting  the  expansions  for  f 
and  w  into  the  stream  function  equation,  one  gets 


n 


o 


Similarly,  because  velocity  scales  become  large  as  Ra  -+■  °°,  the  viscous 
terms  in  the  energy  equation  can  also  be  neglected.  Assuming  that 
temperature  is  of  order  one  in  the  core,  from  the  anticipated  match 
with  the  wall  boundary- 1 ayers ,  Yq  ~  1 .  Thus,  to  first  order,  the 
inviscid  core  expansions  become: 


uc  =  Ra1/2  uQ  +  0(1) 
vc  =  Ra1/2  vQ  +  0(1) 
wc  =  Ra1/2  wQ  +  0(1) 
fc  =  Ra1/2  fQ  +  0(1) 

Tc  *  T0  +  0(Ra'1/2)  .  (4.2) 


With  these  expansions,  the  inviscid  governing  equations,  for  a  finite 


gap  inviscid  core,  take  on  the  following  form  (to  first  order  in  the 
1  imit  as  Ra  -*■<»): 


Energy: 
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St  0  3r 
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Vorticity : 
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Stream  Function: 
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(4.3) 


(4-4) 


(4.5) 
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where  t  =  Ra  '  t. 

With  the  above  time-scaling  factor,  the  local  acceleration  terms 
are  retained  in  the  governing  equations.  Although  steady-state 
analytical  solutions  are  sought,  the  unsteady  terms  enable  one  to  solve 
the  equations  numerically  by  marching  in  a  time-accurate  fashion.  Thus 
the  capability  exists  for  resolving  either  steady  or  unsteady  effects 
within  the  resulting  flow  field.  Also  of  importance,  the  Rayleigh 
number  dependence,  which  was  one  of  the  key  dimensionless  parameters 
in  the  2-D  N-S  equations,  has  been  completely  factored  out  in 
Eqs.  (4.3)  and  (4.4). 


Considering  the  limit  as  G  0,  the  curvilinear  metrics  drop  out 
and  new  scales  are  obtained  in  order  to  retain  all  the  nonlinear 
convective  terms  while  satisfying  the  stream  function  equation: 

u  -  G  u  +  O(G^) 
o  o 

v0  ~  v0  +  0(G) 
w0  ~  G'1  wQ  +  0(1) 
f„  ~  G  fQ  +  0(S2) 

T0  -  T0  +  0(G2)  (4.6 

o 

Neglecting  terms  of  0(G)  and  higher,  relative  to  terms  of  0(G), 

Eqs .  (4.3)  through  (4.5)  become: 


Energy: 


„  3T  ^  3T 
+  uo/  +  vor=° 


(4.7 


Vorticity : 


3w 

Pr  {-4 
3t 


+  U 


o  3r 


+  v 


3w 
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O  3^ 


}  =  s  i 


3r 


(4.8 


Stream  Function: 
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In  these  equations,  the  G  dependency,  which  was  again  important  in  the 
2-D  N-S  equations,  has  also  vanished.  The  stream  function  equation 
appears  to  be  in  boundary-layer  format,  as  should  also  be  the  case  for 
the  energy  and  vorticity  equations  when  the  viscous  terms  are  included. 
The  inviscid  core  equations  can  satisfy  the  flow  tangency  condition, 
but  cannot  satisfy  the  physical  no-slip  condition  that  must  be  met 
at  the  walls,  where  the  viscous  terms  play  their  important  role.  At 
this  point,  one  gathers  from  Eq.  (4.8)  that  the  gravitational  terms 
have  reduced  to  simply  sin<|J  3TQ/3r,  as  G  -*•  0.  This  first-order  term 
is  essential  in  driving  the  buoyancy-induced  flow  field,  and  its  effect 
seems  just  as  significant  in  the  adjacent  boundary- 1 ayers ,  as  will  be 
shown  in  the  next  section.  Note,  however,  that  the  gravitational 
term  proportional  to  cosiji  3Tq/3>Jj  comes  in  as  a  second-order 
effect  within  the  small-gap  inviscid  core  solution.  But  its  effect 
should  prove  negligible  (or  of  much  less  importance)  within  the 
boundary-layers  since  being  consistent  with  boundary-layer  theory, 

3 T / 3 d-'  «  3T/3r . 


4.2.  The  Boundary-Layer  Expansion 
To  incorporate  the  higher  order  viscous  derivatives,  the  inner 
radial  coordinate  is  stretched  as 


r  -  a  r  G  a 


where  5  is  the  inner  boundary -layer  thickness.  Then,  in  its  t ra ns 4 •-•me- 


-•  j.--  '>■ 
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state,  the  tangential  velocity  (v  =  ^  |£)  within  the  boundary-layer, 


Vg^,  becomes  (for  a  finite-gap  size): 


BL  '  5  3R 


(4.11) 


and  due  to  asymptotic  matching  principles,  one  must  have 


VBL  ~  vc  '  0(Ra 


(4.12) 


For  this  condition  to  hold,  the  stream  function  variable  within  the 
boundary-layer  must  scale  as 


fBL  ~  0(5  Ra1/2)  . 


(4.13) 


-1 

The  radial  velocity  (u  =  - y-  ^r-)  then  goes  as  f  (for  finite-gap), 

G(r  +  1)  ^ 

or  G 


ugL  0(5  Ra 


(4.14) 


From  the  definition  of  vorticity,  one  gets 


1  /  2 

V  -HRa  /:) 


(4.15) 


a  s  1 1  i .  t  * *mpe  r  a  t  a  r“  ■  s 
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expansions  are: 
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1/2  F  + 


-  6  Ra,/L  F 


TBL  =  To  +  ••• 


Using  the  above  boundary-layer  expansions  together  with  the  fact  that 
6  0  in  the  limit  as  Ra  -  one  can  reduce  the  2-D  N-S  equations  to 
the  following  form. 


Energy: 
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Stream  Function: 
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and  again,  t  =  Ra  t. 


These  equations  clearly  bring  forth  the  scaling  required  for  the 
boundary- layer  thickness,  that  is 


Ra 


•1/4 


(4.18) 


where  the  1/4  power  is  typical  for  natural  convective  laminar-type  flows. 
Using  Eq.  (4.18),  the  boundary-1 ayer  expansions  become: 


BL 

Ra 

1/4 

Uo  + 

0(Ra 

-1/4 

BL  = 

Ra 

1/2 

Vo  + 

0(1) 

BL 

Ra 

3/4 

Wo  + 

0(  Ra 

!/4) 

BL 

Ra 

1/4 

F  + 
0 

0(Ra 

-1/4 

BL  = 

T 

0 

+  i 

3(Ra" 

1/2) 

(4.19) 


Physically  these  expansions  make  sense,  since  within  the  boundary- 1 ayer , 
vorticity  should  be  much  larger  than  the  stream  function  or  velocity 
components,  and  the  tangential  velocities  should  be  greater  than  the 
radial  velocities  for  this  particular  flow  geometry.  Also,  following 
the  same  aforementioned  steps,  the  outer  boundary-layer  thickness  can 
be  shown  to  scale  in  the  same  manner  as  the  inner  one. 

Now,  as  G  -  0,  the  boundary-layers  must  merge  at  some  point 
giving 
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25  =  G  a  .  (4.20) 

Knowing  this  relation,  G  must  then  scale  as  <5  or 

G  Ra‘1/4  .  (4.21) 

Substituting  Eq.  (4.21)  into  Eq.  (4.6)  and  placing  the  result  into 
Eq.  (4.2),  yields  the  inviscid 

u 

c 

v 

c 

w 

c 

f 

c 

T 

c 

Comparing  Eq.  (4.22)  with  (4.19),  one  can  claim  that  at  the  point 
where  the  boundary- 1 ayers  merge,  the  'outer'  or  core  expansions  do 
indeed  match  the  'inner'  or  boundary- 1 ayer  expansions. 

4 . j .  The  Analytical  Cel  1 -Development  Regime 
Since  the  radius  ratio  b/a  is  of  0(1)  for  narrow  gaps  and  Eqs. 


core  expansions  as 


D  1/4 

=  Ra  u  +  . . . . 
o 


=  Ra1/2  v  ♦  . . . . 
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-  Ra3/4  w  ♦ 
o 


Ra1/4f  ♦ 
o 


•  T  f  . .  .  . 
o 


(4.22) 


(4.17)  are  independent  of  G,  the  boundary- 1 ayer  expansions  do  not  have 
to  be  rescaled  to  avoid  losing  any  nonlinear  physics  as  G  *  0. 


Therefore,  it  can  now  be  inferred  that  the  dependent  variables  in  the 
governing  N-S  equations  will  scale  the  same  as  in  the  boundary-layer 


for  the  double  limit  of  Ra  -*■  °°  and  G  -*■  0.  Thus,  we  take 


and  the  expansions 


G  =  Ra"1/4  G 
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l/2) 

(4.23) 


(4.24) 


With  Eqs.  (4.23)  and  (4.24)  in-hand,  the  2-D  N-S  equations  collapse  into 
Cartesian-like,  boundary-layer  type  equations.  In  final  form,  they 
become : 


Energy : 

~r2  dJ  A  x  r  3f  3T  ,  3f  3T,  _  1  32f 
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Vorticity : 


Pr  ^  +  G 


3t 


pP  r  if  3w  ,  3f  3w  , 
r  l_  3*3r  3r  I  J 


2~ 

=  Pr  +  G  sin  4< 

3r  8r 
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Stream  Function: 
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with  related  boundary  conditions  of: 


T  (0,*)  -  1,  T  (1,*)  =  0 


f  (0,4<)  =  f  (1  ,’!')=  0 


2~ 

w  (O.IP)  =  w  (1,¥)  = 

G^  dr£ 


r=0,l 


(4.25 

(4.25 

(4.25 


where  again,  G,  the  scaled  gap,  is  defined  as 


G  = 


(4.26 


and  the  time-term  is  also  scaled  the  same  as  in  the  boundary-layer,  or 


t  = 


(4.261 
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Equations  (4.25)  represent  the  key  governing  equations  for  arbitrary 
Prandtl  number  in  the  double  limit  of  Ra  -*•  00  and  G  -*  0.  The  numerical 
method  used  to  solve  these  equations  will  be  described  in  Chapter  5. 
These  equations  can  be  further  simplified  under  certain  limiting 
conditions.  Three  limiting  cases  will  be  investigated;  namely,  the 
G  -*•  0  limit,  plus  the  vanishingly  small  and  infinite  Prandtl  number 
1 imits. 

To  leading-order ,  in  the  limit  as  G  •*  0,  Eqs.  (4.25)  reduce  to: 


Energy : 


Vorticity : 


Stream  Function: 


4  =  0 

3r 


(4.27 


Thus,  the  energy  equation  simplifies  to  the  steady-state  1-D 
conduction  equation,  which  integrates  to 


T  *  1  -  r 


(4.28 


for  the  isothermal  boundary  conditions  of  this  problem.  Also,  from 
the  no-slip  condition,  the  vorticity  and  velocity  components  become 
identically  zero.  Hence,  convective  effects  are  discarded,  and  the 
temperature  field  becomes  completely  determined  by  the  linear 
relationship  of  Eq.  (4.28). 

From  the  above  analysis,  one  can  conclude  that  a  Stokes-layer  or 
a  conducti ve/ vi scous  dominated  flow  field  prevails  as  G  -  0.  Therefore, 
the  multicellular  regime  that  is  associated  with  the  narrow  type  gaps 
(at  high  Rayleigh  number)  most  likely  reverts  back  to  two  cells  and 
then,  finally  to  a  pure  conductive  mode  when  the  velocities  tend  to 
zero  as  G  -  0  (or  G  -  0  for  a  specified  Rayleigh  number).  Contrary  to 
this  result,  many  investigators  (Walton,  1980;  Powe  et  al_.  ,  1971;  Liu 
et  aK  ,  1962)  believed  that  as  the  gap  approached  zero,  a  true  Benard 
type  instability  would  evolve,  where  the  typical  critical  Rayleigh 
number  (Ra.  )  of  about  1,700  is  approached.  Instead,  it  appears  that 

D“d 

after  a  certain  point,  the  transitional  Rayleigh  number  increases  with 
decreasing  gap  size  (to  well  above  1,700),  until  eventually,  a  conduction- 
dominated  flow  results  as  the  gap  approaches  zero  (see  Chapter  6). 

4.4.  The  Perturbative  Solution  to  the  Steady-State 
Fini te-Prandtl  Number  Equations 

From  numerical  results  obtained  in  this  study,  and  results 
discussed  in  Powe  et  a(L  (1971)  and  Rao  et  al .  (1985),  it  can  be  seen 
that  the  narrow-gap  solutions  for  flows  prior  to  multicellular 
transition  (and  even  beyond  to  a  certain  extent)  behave  in  a 
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seemingly  steady-state  manner.  Assuming  the  steady-state  condition 
and  using  the  leading-order  solution  of  Eq.  (4.27),  the  convective  terms 
in  the  vorticity  equation  of  (4.25)  become  zero.  From  this,  w  must 
then  scale  as  G  in  order  to  balance  viscous  and  buoyancy  effects,  and 
upon  entering  the  stream  function  equation  of  (4.25),  one  sees  that 
f  G3.  This  logic  can  be  repeated  in  a  successive  manner  (without 
neglecting  cross  terms)  until  the  following  asymptotic  expansions  are 
obtained  for  the  steady-state  case: 


T  =  (1-r)  +  G4T1  +  G8T2  +  G'2T3  +  0(G16) 
w  =  G  w.|  +  G^w^  +  G^w^  +  G  3w^  +  O(G^) 
f  =  G3f1  +  G7f2  +  G1]f3  +  G15f4  +  0(G19) 

while  the  related  boundary  conditions  are: 


(4.29a) 


T(0,<P)  =  1,  TO,*)  =  0 

f(0,*)  =  fO,*)  =  0 

If  (0,*)  =  If  0,*)  =  0  .  (4.29b) 


All  of  the  important  physics  are  included  in  the  first  two  terms  of 
the  expansions  in  Eq.  (4.29a);  that  is,  the  gravitational,  viscous  and 


convective  effects  are  accounted  for. 


Substituting  Eq.  (4.29a)  into  Eq.  (4.25)  and  matching  terms  of 
like  coefficients,  yields  an  infinite  set  of  steady  uncoupled  linear 
partial  differential  equations.  These  equations  take  on  tM  following 
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Vorticity : 
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Stream  Function: 
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Equations  (4.30)  -  (4.38)  were  solved  analytically  with  the  following 
result  (note:  vorticity,  temperature  and  stream  function  were  carried 
out  to  the  first  three  nonzero  terms): 
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where 


A1  =  r2 (r-1  )2 

(4.47) 

B1  =  r(r-l )  +  g- 

(4.48) 

Cl  =  r(2r5  -  6r4  +  5r3  -  1 ) 

(4.49) 

and  A2  through  H4,  which  are  also  functions  of  r  only,  are  given  in 
Appendix  A. 

In  Equations  (4.39)  -  (4.46),  the  stream  function  and  vorticity 
expansions  have  the  same  angular  dependence,  and  to  leading-order, 
are  inversely  proportional  to  Prandtl  number.  For  temperature,  however, 


the  conductive  term  (1  -  r)  and  T-j  are  completely  independent  of  the 
Prandtl  number. 

From  Eq.  (4.39),  it  can  be  seen  that,  to  leading-order ,  the 
stream  function  is  positive  for  all  ^  at  a  given  value  of  r  (for  the  half 
annulus),  and  becomes  zero  at  <l>  =  0,  n.  Therefore,  the  primary  basic 
flow  does  indeed  consist  of  two  counter-rotating  kidney-shaped  cells 
as  discussed  in  the  literature  survey  of  Chapter  2.  It  is  interesting 
to  note  that  the  second  and  third  terms  of  the  stream  function 
expansion  are  proportional  to  sin  2^;  thus,  the  basic  flow  has  the 
potential  of  obtaining  more  than  two  cells  in  the  complete  annulus  - 
especially  when  the  Prandtl  number  is  sufficiently  small  (Mack  and 
Bishop,  1968;  Walton,  1980). 

Before  completing  this  section,  the  asymptotic  expansions  of 
Eq.  (4.29)  will  be  placed  in  a  compact  recursive  notation: 

T  =  (1  -  r)  +  E  G4n  T 
•,  n 

n=l 

00 

p4n-3 

w  =  Z  G  w 
i  n 

n=l 

0° 

f  =  I  G  n_l  f  .  (4.50) 

n  =  l  n 

The  angular  dependent  coefficients  (to  be  denoted  by  F,  W  and  T)  in 
Eqs.  (4.39)  -  (4.46)  can  also  be  written  in  the  following  revealing 


notation : 


sin  V 


-  ;  m=2  to 


F 

m 


OD 


T 


n 


3F 

n 

3<P 


n=l  to  °° 


(4.51  ) 


Since  the  Tn  terms  are  always  functions  of  cos  *•  or  squared  powers 

of  sin  ^  and  cos  the  isotherms  must  remain  symmetric  on  either  half 

of  the  annulus.  Similarly,  since  F  and  W  are  functions  of  sin 

m  m 

sin  or  a  combination  thereof,  the  angular  dependent  coefficients 
will  always  be  equal  in  magnitude,  but  opposite  in  sign  for  each  half 
annulus.  Thus,  due  to  the  above  trigonometric  arguments  (and  since 
symmetry  is  likely  enforced  by  the  concentric  cylinder  geometry), 
symmetric  (bi-cellular  or  multicellular)  solutions  are  implied  for 
either  temperature,  vorticity  or  stream  function.  This  observation  is 
further  supported  by  the  analytical  and  numerical  steady-state 
calculations  obtained  in  this  study  (see  Chapter  6),  along  with  the 
experimental  results  of  Kuehn  and  Goldstein  (1976a)  and  Powe  et  al . 
(1969).  Although  this  analysis  supports  the  result  of  symmetric 
solutions,  especially  at  very  small  gaps,  the  introduction  of  full 
nonlinear  effects  (from  either  the  thermal -energy  or  vorticity  equation 
may  allow  asymmetric  solutions  to  come  into  play.  Therefore,  even 
though  subsequent  numerical  solutions  will  concentrate  on  the  symmetric 
case,  the  asymmetric  possibility  should  not  be  entirely  ruled  out. 


he*  Ter'  jmpare  neat  •'\ir-.,er  predictions  of  the  asymptotic 
met  nop  t,  *^e  r-I'  V:  n  ume r’ ■_  ;  ’  results,  analytical  expressions  for 
we-e  determined.  :'ii'st.  ’nner  and  outer  Nusselt 
njr'Oe1"-.  are  penned  as 


The  first  term  in  Eqs.  (4.54)  and  (4.55),  6”\  represents  the  purely 

conductive  contribution  to  heat  transfer.  The  convective  influence 

~3  ~7 

is  provided  by  the  remaining  two  terms  of  order  G  and  G  . 

In  the  limits  as  Pr  -*■  0  and  Pr  -+  °°  (see  Sections  4.5  and  4.7), 
analytical  expressions  for  the  inner  and  outer  Nusselt  numbers  can 
also  be  obtained.  Hence,  from  forthcoming  equations  (4.76)  and 
(4.106),  one  gets: 

For  Pr  0: 

Nui*  =  Nuq*  =  G-1  (4.56) 


/v 

where  G  is  given  by  Eq.  (4.67),  or 


G  = 


G  Pr 


-1/4 


For  Pr 


*  _  ;-l  ,  r 3  /COS  .  rl  r COS  ,  1  N , 

Nu^  G  +  G  'T440-  G  {  Qfi  (TarTm)} 


96  387701 


+  G'  { 


7  rCOS  r  /  1 


7  rsin2^  /  1 


34560  v 277 


Wt)>  +  G  S  W} 


+  0(G1]) 


(4.57) 


....  *  _  r~  1  r 3  /COS'K  ,  n 7  /COS  2^  ,  -1  u 
Nu^  -  G  -  G  (-|44q)  +  G  {  Q<-  (^^7/lc^^)} 


96  '3374633' 


2  .2 

i  P7  r cos  i  1  \-»  ,  p7  rsin  'P  t  1  \  -i  ,  p / pi  1  \ 

G  34560"  777  ^  G  'Ttt'-'  +  °^G  > 


■34560  171 


(4.58) 


Equation  (4.56)  signifies  the  pure  conduction  result  for  Pr  0,  and 
thus,  does  not  retain  any  convective  influence  to  heat  transfer. 

Equations  (4.57)  and  (4.58)  are  similar  to  those  for  finite  Prandtl 
number,  except  terms  proportional  to  Pr"^  have  vanished. 

4.5.  The  Zero  Prandtl  Number  Limit 
In  this  section,  the  previously  derived  small-gap  (finite-Prandtl 
number)  equations  will  be  used  to  analyze  the  theoretical  effects  as 
Pr  •*  0.  From  a  practical  standpoint,  the  behavior  of  liquid  metals  should 
fall  within  this  low-Prandtl  number  approximation,  and  this  analysis  has 
application  to  the  study  of  liquid  metal-cooled  nuclear  reactors. 

It  appears  from  the  literature  that  Pr  -*■  0  limits  have  not  been  studied 
either  numerically,  analytically  or  experimentally  for  the  horizontal 
narrow-gap  concentric  cylinder  geometry.  Lee  and  Korpela  (1983)  have 
studied  a  similar  type  of  natural  convective  flow  between  vertical 
slots,  where  the  left  wall  is  heated.  They  demonstrated  that  in  the 
narrow  (high  aspect-ratio)  slots,  a  stationary  multicellular  flow 
resembling  an  hour-glass  configuration  is  initiated  near  the  center 
of  the  slot,  and  appears  strongest  for  the  near-zero  Prandtl  number 
cases.  But  the  source  of  this  unique  instability  is  hydrodynamic 
rather  than  thermal,  because  for  Pr  0,  the  large  thermal  diffusitives 
that  result  most  probably  rule  out  any  type  of  thermal  perturbations. 
Realizing  this,  one  can  speculate  that  in  the  vertical  portions  of  a 
narrow  horizontal  annulus,  the  same  type  of  hour-glass  shaped 
hydrodynamic  instability  is  possible.  Furthermore,  Korpela  (1974)  and 


Walton  (1980)  discuss  the  fact  (for  flow  due  to  natural  convection  in 
an  inclined  channel,  see  Chapter  2)  that  for  Pr  >_  .24,  a  true  Benard- 
type  thermal  instability  prevails,  whereas  for  0  <  Pr  <  .24,  a  mixture  of 
thermal  and  hydrodynamic  instabilities  result,  until  finally  at  Pr  =  0, 
the  instability  is  purely  hydrodynamic.  Thus,  one  might  now  look  upon 
the  narrow  horizontal  annulus  geometry  as  a  combination  of  horizontal, 
inclined  and  vertical  channels.  Both  types  of  instability  appear  to 
have  a  definite  role  in  the  cellular  patterns  that  develop  in  this 
geometry,  as  a  function  of  Prandtl  number.  Therefore,  the  small-gap/ 
finite-Prandtl  number  equations  derived  in  Eq.  (4.25)  may  prove  useful 
in  studying  the  cellular-structure  variation  with  Prandtl  number. 

But,  since  the  limiting  equations  for  Pr  0  have  not  yet  been  obtained, 
the  following  analytical  analysis  was  conducted. 

To  begin,  let  Pr  0,  and  assume  again  that  temperature  is  of 
0(1)  from  the  boundary  conditions.  Then,  assume  that 

w  =  A.|  W  +  .... 
f  A2  ^  .... 

and  consider  the  stream  function  equation  in  (4.25)  to  get 

A2  ~  G2  A]  .  (4.61) 

Using  Eqs.  (4.59)  and  (4.60),  the  vorticity  equation  in  (4.25)  becomes: 


(4.59) 

(4.60) 


SPrXl  ^  +  pr  G  Vl  <-  fir  fF  +  fF  I5> 


Pr  X.  +  G  sin  ♦ 

1  3r2  9r 


or 


G2PrX, 


Pr  G  X2X1  -  Pr  X-j  -  G  , 


and  in  order  to  retain  the  viscous  term. 


(4.62) 


(4.63) 


A  ~  — 

A1  Pr  ’ 


then  from  Eq.  (4.61), 


X  -li 

a2  Pr 


Equations  (4.64)  and  (4.65)  are  then  used  in  Eq.  (4.63)  to  obtain 


G  -  0(Pr1/4) 


and  t  ~  G  ,  which  yields  from  Eq.  (4.66a) 


t  ~  0(Pr1/2) 


(4.f^L' 


Thus,  let 


G  =  G  Pr 1 


t  =  t  Pr' 


Hence, 


X1  ~  Pr 


Xg  ~  Pr" 


With  these  results,  the  following  expansions  can  be  written: 


w  ~  Pr_3/4  W  +  0(Pr1/4) 
f  ~  Pr"1/4  F  +  0(Pr3/4) 


T  ~  T  +  0(Pr)  . 


Also,  Eqs.  (4.67a)  and  (4.67b)  can  be  expressed  in  terms  of  the  key 
dimensionless  variables,  namely: 


S.PT-V4J. 

i  -  ^  t  -  (£),/2 1 


These  new  scales,  when  written  in  the  above  form,  clearly  display 
their  unique  dependency  on  Ra,  Pr,  G  and  t.  Upon  substituting  Eqs. 
(4.70),  (4.71)  and  (4.72),  Eqs.  (4.25)  reduce  to  the  following  set 
in  the  limit  as  Pr  -*■  0: 


Energy : 


Vorticity: 


at 


a2T 


=  o 


g2  ^  +  g  {■ 


aF  aw  a£  aw, 
i  3r  3r  ar 


ar  3r 


(4.7: 


(4.7- 


Stream  Function: 


4=  s2  w 

3r 


(4.7! 


The  energy  equation  in  this  case  also  reduces  to  the  steady-state 
conduction  equation.  For  the  boundary  conditions  considered,  the 
linear  temperature  distribution  across  the  gap  takes  on  the  following 
simple  solution: 


T  =  1  -  r  . 


(4.7i 


Thus,  as  Pr  -*■  0,  only  concentric  isotherms  should  result,  regardless 
of  the  flow  field  that  develops.  The  energy  equation  has  now 
completely  uncoupled  itself  from  the  vorticity  equation,  and  since 


Eq.  (4.76)  gives  3T/3r  =  -1,  the  final  two  governing  equations,  which 
are  directly  independent  of  Prandtl  number,  become: 


Vorticity: 


^  +  G  {- 

at 


3F  3W  3F  aw,  _  a2w 
W  ar  ar  ar  "  77 


6  sin 


(4.77) 


Stream  Function: 


a2F  _  :2 
-  G 

3r 


(4.78) 


with  boundary  conditions 

F(0,>|J)  =  F(l,40  =  0  (4.79) 

|f  (0,*)  =  |p  (1»  =  0  .  (4.80) 

These  equations  will  be  solved  numerically  in  order  to  determine  if  a 
hydrodynamic-type  instability  (similar  to  that  previously  described)  does 
indeed  exist  for  this  flow  geometry  (see  Chapter  6  for  a  discussion  of 
these  results). 

4.6.  The  Perturbation  Solution  to  the  Steady-State 
Zero  Prandtl  Number  Equations 

Assuming  steady-state  and  upon  eliminating  W,  one  gets 


kw. 


93 


/\ 

G 


r  3F  83F  ,  3F  33F  , 


sin  Q 


(4.82) 


The  logic  used  to  derive  the  expansions  in  Eq.  (4.29)  when  G  0  can 
again  be  used  to  derive  the  expansions  for  the  G  -*  0  limit.  Thus, 
considering  Eq.  (4.82),  the  stream  function  expansion  should  behave  as 


/SO  A1  ^11  C 

F  =  GJF1  +  G7F2  +  G  F3  +  0(G  j 


(4.83) 


and  from  Eq.  (4.78),  the  vorticity  expansion  becomes 

W  =  G  W-|  +  G5W2  +  G9W3  +  0(G13)  .  (4.84) 

Notice  that  the  power  of  the  coefficients  in  the  above  expansions  is 
identical  with  those  in  Eqs.  (4.29),  as  should  be  the  case.  However, 
in  deriving  the  asymptotic  expansions  of  Eqs.  (4.29),  the  buoyancy 

~  OX 

term,  G  sin  V  (in  Eq.  4.25b),  maintained  the  same  order  as  the 
viscous  term  for  all  orders  of  the  vorticity  expansion.  But  for  the 
Pr  -+■  0  limit,  the  buoyancy  term  comes  directly  into  play  only  in  the 
leading  order  W-j  term,  where  the  viscous  term  is  balanced  by  the 
buoyancy  force;  while  in  the  other  terms  of  the  vorticity  expansion, 
the  convective  forces  are  equated  with  the  viscous  ones  and  the 
buoyancy  effects  contribute  only  indirectly  through  F^  and  . 

Substituting  expansion  (4.83)  into  Eq.  (4.82),  and  matching  terms 
of  like  coefficients,  yields  the  following  steady-state  result: 


3F,  d  F,  3F,  a  F 


Integrating  Eq.  (4.85),  one  gets 


4  3  2 

F,  =  Jr  sin  *  +  c,  +  c2  +  c3r  *  c4  , 


(4.87 


while  the  constants  of  integration  are  obtained  from  the  boundary 
conditions  of  Eqs.  (4.79)  and  (4.80)  to  give: 


F,  ■  r2  (r-1)2 


(4.88 


Differentiating  Eq.  (4.88)  twice  with  respect  to  r  gives 


(4.89 


Using  Eq.  (4.88),  Eq.  (4.86)  becomes 


—  =  {6r5  -  15  r4  +  14  r3  -  6  r2  +  r}  2-s~"  24>  .  (4.90 

_  T  O 


Integrating  Eq.  (4.90)  and  applying  boundary  conditions  (Eq.  4.79)  and 
(Eq.  4.80)  yields 


KH*1 
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and  follows  as 

W2  =  •§i^  (D2)  (4.92) 

where  B2  and  D2  are  the  same  as  those  in  Eqs.  (4.42)  and  (4.43),  and 
again,  are  given  in  Appendix  A. 

The  equations  for  and  Wg  take  on  forms  similar  to  those 

in  Eqs.  (4.39)  through  (4.43),  except  that  in  the  former,  the  Prandtl 
number  dependency  and  the  coefficients  A2  and  C2  completely  vanish  for 
the  limit  of  Pr  -*■  0.  Then,  by  neglecting  terms  A3,  C3,  E3  in  Eq.  (4.45) 
and  A4,  C4,  E4  in  Eq.  (4.46)  for  the  limit  as  Pr  -*■  0,  one  easily 
obtains  solutions  for  F^  and  W3,  respectively. 

Since  F2,  F^  and  W2,  W3  involve  terms  proportional  to  sin  2^,  it 
can  be  reasoned  that  these  terms  must  contribute  to  the  potential 
multicellular  transition  of  the  flow  field  occurring  at  some  critical 

A 

value  of  G.  By  performing  a  simple  order-of-magnitude  analysis,  one 
can  estimate  when  these  terms  become  of  equal  importance  with  the 
leading  order  expansion  terms.  From  this  result,  the  onset  of 
multicellular  flow  can  be  approximately  predicted.  Carrying  out  this 
analysis,  one  finds  that  the  sum  of  G^F^  and  G^F^  achieves  the  same 
order  as  G^F^ ,  when  G  =  7.2.  Similarly,  for  vorticity,  G^W^  and  G^W^ 
achieve  the  same  order  as  GW-]  at  G  =  7.5.  Hence,  the  critical  value 
of  G  appears  to  be  about  7. 
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For  the  pretransi tional  cases,  the  first-order  expansions  of 
Eqs.  (4.83)  and  (4.84)  provide  good  approximations  in  describing  the 
behavior  of  stream  function  and  vorticity  (respectively)  for  the  limit 
as  Pr  ■+  0.  Since  it  is  likely  that  terms  beyond  and  W3  play 
dominating  roles  in  triggering  the  instability,  the  critical  value  of 
G  might  be  lower  than  that  predicted  by  the  above  order-of-magnitude 
analysis.  To  properly  estimate  this  value,  the  full  effects  of 
nonlinearity  must  be  taken  into  account.  This  will  be  done  in  Chapter  6 
when  the  equations  are  solved  numerically.  Considering  the  above 
discussion,  one  obtains  the  following  relations: 

F  =  gV  (r-1)2^  (4.93) 

W  =  G  {r  (r-1 )  +  1}  .  (4.94) 


Equations  (4.93)  and  (4.94)  represent  analytical  profiles  for  the 
steady-state  pretransi tional  stream  function  and  vorticity,  which 
are  valid  throughout  the  annulus  (0  -  2m).  Comparisons  to  related 
numerical  results  will  be  given  in  Chapter  6.  From  these  relations, 
it  can  be  seen  (for  the  half-annulus)  that  the  vorticity  given  by 
Eq.  (4.94)  takes  on  positive  values  within  the  boundary-layers  (for 
r  near  zero  and  one),  and  negative  values  within  the  inviscid  core. 

The  stream  function  (from  Eq.  4.93)  always  remains  non-negative,  since 
the  r  terms  are  squared.  This  reasoning  can  be  extended  throughout  the 
entire  flow  field,  at  least  for  the  bicellular  case. 


97 


Referring  to  Eqs.  (4.93)  and  (4.94),  it  can  be  seen  that  any 
resulting  bicellular  profiles  must  be  symmetric  with  respect  to  both 
the  vertical  and  horizontal  centerlines,  or  at  ^=0,  j,  it  and  .  This 
suggests  that  for  the  Pr  0  limit,  the  pretransitional  flow  fields 
are  completely  symmetric,  although  the  possibility  exists  for  the 
occurrence  of  asymmetric  solutions  due  to  nonlinearity,  especially 
after  the  onset  of  instability.  In  addition,  for  the  half-annulus 
basic  cell -structure,  the  stream  function  achieves  its  maximum,  and 
vorticity  its  maximum  or  minimum  (depending  on  r),  at  =  j  for  any 
particular  r  location. 

Furthermore,  because  the  extremum  stream  function  location 

occurs  at  (or  very  near)  =  j,  for  Pr  -*■  0,  and  since  the  buoyancy 
3T 

term  sir,  ^  reduces  to  -sin  V  for  this  case,  one  can  speculate  that 
d  r 

the  upward  or  downward  shift  in  the  stream  function's  position  (caused 

by  an  increase  in  Rayleigh  number)  is  probably  dependent  upon  the 
3T 

magnitude  of  and  the  type  of  fluid  used  -  or  the  value  of  Prandtl 

d  r 

number. 


4.7.  The  Infinite  Prandtl  Number  Limit 
To  properly  complete  the  analytical  analysis  of  this  particular 
flow  geometry,  the  other  end  of  the  Prandtl  number  spectrum  will  be 
considered.  Thus,  simplified  governing  equations  are  sought  for  the 
Pr  -o  limit.  Many  important  fluids  (such  as  glycerin  and  most  oils) 
possess  a  large  Prandtl  number  and  should  be  fairly  accurately  described 
by  the  infinite  Prandtl  number  approximation. 
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Let  Pr  -*■  00  and  assume  again  that  temperature  is  of  0(1);  then  let 


w  =  n-j  w  +  —  • 

f  =  ^  +  •  •  •  • 

T  =  T  +  ....  (4.95 


Substituting  the  above  into  the  finite-Prandtl  number  equations  of 
(4.25)  gives 


Energy: 


2; 


r2  nv.  01"  .  Dv.  :r  9f  9T  .  9f  9Tn  _  9  T  ,,  Qf- 

G  .  Pr  +  ^  .  Pr  G[-  w  ^  +  -  w]  -  (4.96 

9t  9r 


Vorticity : 


o  9T  D  pr  9f  9w  9f  9w-| 

Pr  nl  ^  +  n2nl  '  Pr  G[‘  S  IF  +  sF 


2- 

3  w  x  ,h  3T 


_-nipr^+Gsin^ 


(4.97 


Stream  Function: 

2- 

ru  — jr  =  ru  0^  w  .  (4.98 

4  9r 

In  the  limit  as  Pr  -*•  »,  for  6  and  T  of  0(1),  in  the  viscous  term 
of  the  vorticity  equation  must  go  as  ^  so  as  to  maintain  the  same 
order  as  the  buoyancy  term,  which  is  also  regarded  as  0(1).  Now  the 


o.y.v.v. 


x*->v.y  v- 
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stream  function  equation  must  be  satisfied;  thus  from  Eq.  (4.98),  one 
obtains 


n2~  nl 


1_ 

Pr 


(4.99) 


and  this  expression  for  also  satisfies  the  energy  equation  of  (4.96) 
at  steady-state.  However,  the  unsteady  time  term  in  the  energy  equation 
can  be  retained  by  scaling  t  as 


t  =  Pr  t 


(4.100) 


or  from  Eq.  (4.26b)  one  gets 

-  1  -  Ra1/2 

(5pr-)t  •  (4.101) 

Hence,  for  the  Pr  -*■  °°  limit,  the  following  expansions  result: 

_1  -  _2 
w  =  Pr  w  +  0(Pr  ) 

f  =  Pr"1  f  +  0(Pr"2) 

T  =  T  +  0(Pr-1 )  .  (4.102) 


With  Eqs  (4.100)  and  (4.102)  in-hand,  the  governing  equations  become 
(in  the  limit  as  Pr  -*■  <*>): 


Energy: 


(4.103) 


(4.104) 


(4.105a) 


T(0,40  =  1,  T(l»  »  0 
f(0,<P)  =  f(l,*)  =  0 

and  fj  (0,*)  =  §£  (1,V»)  =  0  .  (4.105b) 

By  letting  Pr  -*■  <»,  the  nonlinear  convective  terms  in  the  vorticity 
equation  have  completely  dropped  out.  However,  since  the  nonlinear 
terms  in  the  energy  equation  still  remain,  the  decoupling  effect 
between  the  energy  and  vorticity  equations  is  no  longer  possible,  as 
it  was  with  the  Pr  -*■  0  limiting  solution.  Also  in  the  above  equations, 
the  Prandtl  number  dependency  has  again  been  eliminated. 


Similar  to  the  finite-Prandtl  number  equations,  the  G  -*■  0  limit 
yields  a  purely  conductive-dominated  flow  field,  where  again,  the 
temperature  profile  reduces  to  the  1-D  steady-state  conduction  case. 
Hence,  steady-state  expansions  can  once  again  be  found  by  proceeding 
with  the  G  -*■  0  solution.  As  expected,  the  expansions  that  result  take 
on  the  same  form  as  those  in  Eq.  (4.29a): 

T  =  (1-r)  +  G4  T1  +  G8  T2  +  0(G12) 
w  =  G  W|  +  G8  v>2  +  0(G8) 

f  =  G3  f1  +  G7  ?2  +  0(G1])  .  (4.106) 

Perturbative  solutions  to  the  infinite  Prandtl  number  equations  are 
readily  found  by  neglecting  the  Prandtl  number  dependency  and  the 
coefficients  of  O(^)  relative  to  those  of  0(1),  in  Eqs.  (4.39)  through 
(4.46).  This  approach  is  similar  to  that  used  in  the  Pr  -►  0  limit, 
but  with  opposite  terms  neglected. 

Finally,  then,  in  the  limits  of  small  and  very  large  Prandtl 
number,  the  governing  equations  for  the  narrow-gap  limit  have  become 
much  simpler.  In  the  former,  the  nonlinear  terms  in  the  vorticity 
equation  (which  could  give  rise  to  hydrodynamic  instabilities)  dominate, 
whereas  in  the  latter  case,  the  nonlinear  terms  in  the  energy  equation 
(which  should  give  rise  to  thermal-type  instabilities)  dominate.  Thus, 
by  examining  these  two  limiting  conditions,  the  crucial  mechanism  for 
causing  potential  natural  convective  instabilities  between  narrow 


concentric  horizontal  isothermal  cylinders  is  clearly  set  forth. 

Simplified  governing  equations  for  Pr  -►  0  and  Pr  -*■  <=°  can  also  be 
obtained  from  the  2-0  N-S  equations,  which  include  both  buoyancy 
terms  and  are  valid  for  arbitrary  gap  size  and  Rayleigh  number.  These 
limiting  conditions  are  explored  and  discussed  in  Appendix  B. 

Although  the  material  in  this  chapter  is  very  analytical  in  nature, 
many  of  the  equations  derived  (with  the  exception  of  those  relating  to 
the  infinite  Prandtl  number  limit)  will  also  be  solved  numerically  in 
Chapter  6.  Hence,  useful  comparisons  to  the  analytical  results  can 
be  made,  and  the  critical  effects  of  the  nonlinear  terms  can  be 
interpreted. 


103 


5.  NUMERICAL  ANALYSIS 

This  chapter  is  divided  into  two  sections.  The  first  describes 
the  finite-differencing  method  used  in  solving  the  coupled  set  of 
governing  partial  differential  equations.  Two  systems  of  equations 
will  be  considered:  the  2-D  unsteady  elliptic  Navier-Stokes  equations 
derived  in  Chapter  3,  and  the  2-D  unsteady  boundary-layer  equations 
developed  in  Chapter  4.  The  second  section  deals  mainly  with 
describing  the  numerical  procedure  employed  in  determining  the  onset 
of  an  instability,  either  thermal  or  hydrodynamic,  for  both  the  2-D 
Navier-Stokes  and  boundary-layer  equations.  Procedures  for  capturing 
related  hysteresis  behavior  will  also  be  discussed.  In  addition, 
important  computational  details  will  be  provided  at  the  end  of  this 
section. 

The  essential  features  of  the  numerical  method  were  adopted  from 
Prusa  (1983).  However,  several  important  modifications  were 
incorporated  into  his  scheme  in  order  to  enhance  numerical  stability 
and  efficiency  when  solving  the  system  of  equations  in  the  high  Rayleigh 
number/small -gap  multicellular  flow  regime.  Of  particular  importance 
was  the  handling  of  the  nonlinear  convective  terms,  in  which  two  types 
of  representations  were  employed:  a  corrected  second-order  central 
difference  scheme  for  the  2-D  Navier-Stokes  equations,  and  a 
corrected  second-order  upwind  scheme  for  the  boundary-layer-like 
equations.  Also,  the  unsteady  form  of  the  equations  was  used  to  take 
advantage  of  achieving  steady-state  in  a  more  stable  manner  and 
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allowing  the  opportunity  to  capture  unsteady  behavior  when  solved  in 
a  time-accurate  fashion. 

5.1.  The  Numerical  Method  of  Solution 

The  2-D  Navier-Stokes  equations  outlined  in  Eq.  (3.1)  are 
discretized  using  finite-differencing  techniques.  This  converts  the 
equations  into  an  algebraic  system  that  is  well-suited  for  computation 
on  a  high-speed  computer.  In  the  transformed  plane,  the  computational 
domain  was  divided  into  a  number  of  cells.  This  cellular  mesh  was 
formed  by  the  intersection  of  a  set  of  radial  lines  with  a  set  of 
circular  arcs  concentric  with  the  boundary,  r  =  1.  The  grid  nodes  are 
variably  spaced  and  are  defined  by  the  intersection  of  these  radial 
lines  and  concentric  arcs.  The  variable  increment  nodes  are  depicted 
in  Figure  5.1  by  using  an  extended  version  of  Southwell's  (1946) 
notation  in  the  immediate  neighborhood  of  a  typical  node.  Initially, 
the  complete  annular  flow  field  (0  to  2ir)  was  resolved  numerically. 
However,  if  symmetric  cellular  patterns  resulted,  vertical  symmetry  was 
assumed  and  only  the  half-annulus  (0  -  tt)  was  computed  for 
subsequent  calculations. 

The  use  of  variable  increments  allowed  for  the  concentration  of 
grid  nodes  in  areas  of  large  gradients  within  the  computational 
domain,  such  as  the  inner  and  outer  boundary-1  ayc-rs  and  the  thermal 
plume  region  near  the  top  of  the  inner  cylinder,  -lowever,  it  should 
be  noted  that  the  use  of  variable  increments  mav  result  in  a  loss  of 
formal  truncation  error  (Roache,  1976),  where  the  loss  becomes  less 


106 


severe  in  the  finer  portions  of  the  mesh  and  more  severe  in  the  coarser 
areas  (this  result  is  discussed  further  in  subsection  5.1.2).  A 
convenient  2-D  coordinate  variable  mesh  routine  was  developed  in 
Prusa  (1983),  in  which  the  variable  nodal  positions  were  computed  by 
using  a  smooth  fourth-order  polynomial  stretch  that  imposed  a  zero 
gradient  increment  variation  at  the  edges  of  the  mesh.  This  routine, 
with  minor  changes  to  accommodate  the  complete  (0  -  2tt)  annulus,  was 
used  in  this  study  to  generate  the  various  computational  meshes.  The 
program  was  relatively  simple  to  use  and  proved  quite  valuable  in 
resolving  the  multicellular  flow  field. 

5.1.1.  Variable  increment  finite-difference  formulas 

The  formal  Taylor  series  approach  will  be  used  to  develop  variable 
increment  finite-difference  approximations  of  the  derivatives  in  the 
governing  equations.  The  dependent  variables  T,  w  and  f  will  be 
denoted  by  <p ;  k^  and  denote  the  increments  in  the  directions  of 
increasing  and  decreasing  angular  variation,  from  a  given  discrete 
point  ( r .j ,  'J'j)  (see  Figure  5.1).  Using  Southwell’s  notation,  the 
Taylor  series  expansions  of  <t>2  and  <j>^  (r^,  4^)  are: 


a,  =  *  +  k  +  lljfc  k  2  +  k3  + 

*2  *  0  f  2  3*2  0  f  6  3*3  0  f 


4  =  *  -  M  k  +  k  2  _  Ilk  k3  + 

*4  *0  3<P  q  u  2  3^2  Q  kb  6  3lj;3  Q  kb 


(5.1a) 

(5.1b) 


Equation  (5.1a)  can  now  be  solved  explicitly  to  give  a  forward- 


difference  representation  of  3<J>q/ : 


(5.1c) 


while  Eq.  (5.1b)  yields  the  backward-difference  form: 


+  0(kb) 


(5. Id) 


Subtraction  of  Eq.  (5.1b)  from  Eq.  (5.1a)  gives  the  (second-order-like) 
centered-difference  result: 


dip 


0 


~  p  p 

r-rr +  0(kf "  kb)  +  0(kf  >  kfkb’  V  • 


( 5 . 1  e ) 


A  second  linear  combination  of  Eqs.  (5.1a)  and  (5.1b)  gives  the 
centered  second  derivative: 


d2<p 


0 


kb^2  ”  ^kf  +  kb^0  +  ^f^4 

2  (kf  + 


+  0(kf  -  kb)  +  0(k2 ,  kfkb,  k2) 


(5. If) 


The  counterpart  radial  expressions  for  Eqs.  (5.1c)  through  (5. If) 
are  easily  obtained  by  substituting  r  for  4),  and  4>-| ,  4>^  anc^  hb  for 
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<(>2 »  $4  and  kf,  k^,  respectively. 

Finite-difference  formulas  for  the  unsteady  terms  are  found  by 
using  the  standard,  but  generally  stable,  forward-difference  molecule: 


9(j) 

at 


o 


At 


+  o  ( At ) 


(5.2) 


where  n  designates  the  time  level  at  which  the  dependent  variable  <j>  is 
evaluated.  When  marching  in  time  to  steady-state,  the  first-order 
truncation  error  should  not  influence  the  final  result;  hence,  time 
steps  as  large  as  possible  should  be  considered.  However,  when 
unsteady  behavior  had  to  be  resolved  (as  with  the  Pr  -*•  0  case  in  this 
study),  very  small  time  steps  should  be  used  to  avoid  large  truncation 
error,  thus  obtaining  a  time-accurate  solution. 


5.1.2.  Finite-difference  equations  for  the  dependent  variables 

All  spatial  derivatives  in  the  governing  equations  are  second-order 
centrally  (or  upwind)  differenced.  This  includes  the  convective  terms 
which  are  represented  by  a  first-order  upwind  expression  together  with 
a  correction  term  for  second-order-li ke  accuracy.  The  convective  terms 
are  split  in  order  to  maintain  diagonal  dominance  of  the  coefficient 
matrix  at  high  Rayleigh  numbers  (as  with  pure  upwind  differencing), 
thus  providing  stable  convergence  toward  the  desired  centrally 
differenced  (or  second-order  upwind  differenced)  representation. 

5. 1.2.1.  2-D  Navier-Stokes  equations  The  nonlinear  convective 

terms  in  the  2-D  Navier-Stokes  equations  can  be  expressed  as: 


;  S-  r.V.  /. 


.v 


4>n  -  4>o  +  -  (1  +  H)<j>n 

2x|i  -  (X  -  |X|)  {-V1*  J - TTT - 

dr  q  %  nb  nf 


+  (*  +  IM  )  t - h 


*1  "  *0  H1  +  *3  “  {1  +  iH 


hb  +  hf 


(5.3a) 


(5.3b) 


?  3*  _  .  » 0  ~  *4  ,  h  +  k(^4  ~  (1  +  k)^C 
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where  <j>  represents  the  dependent  variables  T,  w,  or  f,  and 


hf  kf 

H  h  »  k  =  IT  * 

nb  Kb 


The  first  terms  within  the  brackets  are  the  first-order  upwind- 
difference  components  and  the  second  terms  are  the  added  corrections 
which  bring  the  differencing  up  to  second-order  accuracy.  The  terms 
within  the  parentheses  (outside  of  the  brackets)  ensure  that  stable 
differencing  ' into-the-wind'  is  enforced. 

Using  the  above  relations,  Eqs.  (3.1a,  b,  or  c)  can  be  manipulated 
into  the  following  general  form: 


A2  {-  AG2  +  ^4  +  (r  +  1)  2  +  2A  +  2u  +  S  =  0  (5.4) 


2A  =  (r  +  £)  (A2  +  AT 


(5.5 


and 
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Pr,  <J>  =  T 
Pr,  <j>  =  w 
0,  4>  =  f 


A2 


1,  =  T 

Pr,  4>  =  w 
1,  4>  =  f 


while 


A 


Pr,  <p  =  T 
1,  4)  =  w 

0,  <J>  *  f 


(5.5 


The  derivatives  in  Eqs.  (5.5a)  and  (5.5b)  are  represented  by  standard 
central -difference  expressions,  as  in  Eq.  (5.1e). 

Employing  Eqs.  (5.3a  and  5.3b)  and  substituting  the  appropriate 
differences  of  Eqs.  (5.1c-f),  results  in  a  general  finite-difference 
equation  for  <j>  =  T,  w,  or  f: 


CrVt>n 


n+1  (A2)  A  G* 


+  C 


n+1 


+  C 


o4>o 


n+1 


Note  that  a  negative  contribution  to  CQ  is  possible  if  hf/hb  or  kf/k 
departs  too  far  from  1.  Thus,  unstable  calculations  may  result  if 
diagonal  dominance  is  lost  as  C  gets  smaller  in  magnitude. 
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’  1,  4>  =  T 
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( 5 . 6  i 


The  derivatives  —  and  §  in  Eq.  (5.6h)  are  evaluated  by  using 

ar  o  w  o 

standard  central -di fference  formulas,  t  denotes  the  time  increment, 
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while  the  superscripts,  n,  n+1 ,  are  as  discussed  in  Eq.  (5.2). 

For  a  fully-implicit  method,  the  correction  term  E,  should  be 
evaluated  at  the  new-time  level,  En+^ .  However,  as  seen  in  Eq.  (5.6g), 
the  correction  term  for  the  2-0  Navier-Stokes  equations  was  evaluated 
at  the  old-time  level,  En.  Although  the  En  method  introduced  a  time 
truncation  error  of  0(At),  it  still  remained  consistent  with  the  rest 
of  the  formulation.  Also,  since  En  uses  known  information  from  the 
previous  time  level,  its  contribution  to  the  off-diagonal  terms  (i^, 

$2,  <J>2  and  was  kept  constant  when  evaluating  flow  properties  at 
the  new-time  level.  This  factor  apparently  helped  to  enhance 
numerical  stability,  which  in  turn,  improved  the  convergence  rate  and 
accuracy  of  the  numerical  method  when  compared  to  the  En+1 
formulation.  This  was  especially  true  for  the  highly  convective  flows 
associated  with  multicellular  development  (see  Appendix  C).  In  comparing 
to  En+^ ,  a  truncation  error  analysis  demonstrated  that  the  En  method 
preserved  the  formal  accuracy  of  the  system  of  equations.  A 
description  of  this  analysis  and  a  tabular  comparison  of  the  two 
approaches  is  given  in  Appendix  C. 

5. 1.2. 2.  Fini te-Prandtl  number  boundary-layer  equations  The 
boundary-layer  equations  of  (4.25)  are  numerically  treated  in  a  manner 
similar  to  that  used  on  the  2-D  Navier-Stokes  equations.  But  since 
Eqs.  (4.25a-c)  are  only  first-order  in  v1,  a  slightly  modified  version  of 

Eq.  (5.4)  results: 

2 

A2  {-  A  G2  H  +  2A  !£  +  2U  +  S  =  0  (5.7) 
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The  first-order  derivatives  in  the  2\i  and  2X  coefficients  are 
approximated  by  central -differences.  Also, 
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Except  for  the  definition  of  A,  the  convective  term  2A  is 
represented  as  in  Eq.  (5.3a).  However,  because  of  the  parabolic  nature 
of  these  boundary-layer  equations,  which  are  first-order  in  4*,  a  central 
difference  representation  of  the  streamwise  convective  terms  is  no 
longer  valid.  Instead,  a  stable  corrected  second-order  upwind- 
differenced  scheme  is  employed.  Namely, 
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(see  Figure  5.1  for  related  nodal  positions).  The  second  term  within 


the  brackets  represents  the  added  correction  to  the  first-order 


differenced  component  needed  to  achieve  a  second-order  upwind- 


differenced  expression. 


Using  the  above  results,  the  boundary-layer  equations  of  (4.25) 
can  be  cast  into  a  finite-difference  form  similar  to  that  of  Eq.  (5.6a): 
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and  lastly, 
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centrally-di fferenced. 


The  first  two  terms  of  Eq.  (5.10g)  represent  the  second-order 
corrections  corresponding  to  central-differencing  in  the  transverse  or 
radial  direction,  whereas  the  last  two  terms  signify  the  corrections 
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necessary  to  obtain  second-order  upwind  differencing  in  the  streamwise 
or  angular  direction.  In  contrast  to  the  2-D  Navier-Stokes  equations, 
the  correction  terms  for  the  boundary-layer  equations  were  evaluated  at 
the  new-time  level,  En+1 .  Since  it  is  usually  the  central -differenced 
streamwise  convective  terms  that  cause  stability  problems  at  high 
Rayleigh  numbers,  the  En  approach  was  adopted  to  alleviate  this 
predicament  for  the  Navier-Stokes  equations.  However,  since  second-order 
upwind  differencing  is  generally  more  stable  than  central -di fferenci ng , 
the  En+^  method  worked  fine  and  remained  stable  for  the  boundary-layer 
calculations.  Note  that  the  En+^  method  also  converged  for  the  Navier- 
Stokes  calculations,  but  the  En  formulation  proved  to  be  more  efficient 
(see  Appendix  C).  Based  on  the  literature  review  and  work  in  this  study, 
it  was  found  that  although  first-order  upwind  differencing  is  extremely 
stable  at  high  Rayleigh  numbers,  at  least  second-order  accuracy  is 
needed  to  resolve  the  multicells. 

5. 1.2. 3.  Zero-Prandtl  number  boundary-layer  equations  In  the 
Pr  -*■  0  limit,  the  energy  equation  reduced  to  simply  T  =  1  -  r,  therefore 
decoupling  itself  from  the  vorticity  equation.  Hence,  only  the  vorticity 
and  stream  function  equations  had  to  be  solved  numerically  in  a  coupled 
manner.  The  following  one-equation  format,  similar  to  Eq.  (5.7),  can 
again  be  used  to  represent  these  two  coupled  equations  as  derived  in 
Chapter  4  (Eqs.  4.77  and  4.78): 


A2  {-  AG2  4  +  ^-4}  +  2>  4  +  2u  M  +  s  =  0 
3t  5r  jr  J 


where 


The  finite-difference  form  of  Eq.  (5.11)  is  identical  to  Eq. 
(5.10a)  for  the  finite-Prandtl  number  boundary-layer  equations,  except 
that  G  is  now  replaced  by  G.  The  same  applies  to  the  coefficient 
equations  for  Cq,  ,  C2,  and  in  Eqs.  (5.10b)  -  (5.1  Of ) ,  but 
with  X  and  y  being  defined  by  Eqs.  (5.12a)  and  (5.12b),  respectively. 
The  En+^  correction  term  is  similarly  represented  by  Eq.  (5.10g), 
while  again  making  the  appropriate  substitutions  for  X  and  y.  Also, 
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5.1.3.  Boundary  conditions 

5. 1.3.1.  2-D  Navier-Stokes  equations  Considering  the  complete 
annulus,  the  finite-difference  form  of  the  boundary  conditions  for  the 
2-D  Navier-Stokes  equations  is: 


Energy : 


(5.14a) 


(5.14b) 


Vorticity : 


w 


n+1 

1J 


2f 


n+1 

ill 


(Gh^  y 


n+1 

WNR,j 


2f 


n+1 

NR-1  ,j 


(GhNR_l ) 


1 


Stream  Function: 


0.0 


(5.14c) 


(5.1 4d ) 


(5.14e) 


0.0 


(5 . 14f ) 


Because  of  continuity  at  zero  and  2tt  radians,  a  computational 
continuous  condition  must  be  defined  such  that: 


rn+l 


A 

i  ,NS 


_  kNS-2Ti,NS  ~  (kNS-l  T  *NS-2;ii,NS-l  T  *NS-l'i,NS-2 


kNS-2^kNS-1  +  kNS-2  ^ 


and  for  Eqs.  (5.14)  and  (5.15),  =  r^  -  r. ,  k^  =  Also, 

the  i  index  refers  to  the  radial  position  with  i  =  1  corresponding  to 
the  inner  cylinder,  r  =  0;  and  i  =  NR  corresponding  to  the  outer 
cylinder,  r  =  1 .  The  j  index  signifies  the  angular  position  with 
j  =  1  referring  to  ^  =  0°,  j  =  NS  referring  to  ^  =  180°,  and  j  =  NS1 
referring  to  =  360°. 

The  symmetry  conditions  for  temperature  were  obtained  from  Taylor- 
series  expansions  about  j  =  1  and  NS,  namely  (for  ^  =  0): 


T  =  j  +11 

i  ,2  1,1  d<P 


k  +iiii 

i  ,1  1  2  d*P2 


1,1 


and  upon  assuming  symmetry,  one  gets: 


II 

dip 


*  o 


i  ,1 


Then,  solving  for  T.  ,  and  substituting  the  second-order  difference 

2*  *  » 

T 

expression  for  (evaluated  at  j  =  1),  the  final  form  of 

Eq.  (5.15a)  is  obtained.  Equation  (5.1 5d )  can  also  be  approximated  in 

a  similar  manner. 


The  formal  Taylor-series  approach  is  also  used  to  derive  the 
vortici ty-wal 1  boundary  conditions.  At  the  inner  wall,  r  =  0,  let: 
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1  ,j 

0,  at  wall 


3f 

3r'i,j 
0,  no-slip 


l  a2f 

+  2  ~ ? 
3r 


U 


0(h1)' 


Then, 


a2f 

a7 


U 


and  from  the  definition  of  vorticity,  one  obtains: 


(Gh, ) 


2  • 


Likewise,  the  condition  for  the  outer  wall  can  be  derived.  According 
to  Roache  (1976),  this  first-order  difference  form  of  the  vorticity 
boundary  condition  is  the  safest  and  most  stable,  and  in  fact,  is 
sometimes  more  accurate  than  other  second-order  approximations. 

5. 1.3. 2,  Finite-Prandtl  and  zero-Prandtl  number  equations  Th 
finite-difference  expressions  for  the  boundary  conditions  of  the 
finite-Prandtl  and  zero-Prandtl  number  equations  are  similar  to  those 
of  the  2-D  Navier-Stokes  equations.  These  expressions  can  be  found 
by  simply  replacing  6  with  G,  G,  and  f,  w,  T,  t  with  f,  w,  T,  t  and 
F,  W,  T,  t,  in  Eqs.  (5.14)  and  (5.15),  respectively.  Note  again  that 
for  the  Pr  -*■  0  case,  the  energy  equation  reduced  to  T  =  1  -  r,  thus 
decoupling  itself  from  the  vorticity  equation.  Therefore,  if  vertical 


symmetry  is  assumed,  the  ^  =  0  constraint  drops  out  of  the  boundary 
conditions  at  ^  =  0  and  <P  =  it. 

Finite-difference  expressions  for  Nusselt  number  and  shear  stress 
are  derived  and  explained  in  Appendix  D. 

5.2.  The  Computational  Procedure 

The  systems  of  coupled  finite-difference  equations  (as  described 
in  the  preceding  section)  were  solved  implicitly  in  time  using  a 
point  iterative  Gauss-Seidel  method  with  underrelaxation.  The 
dependent  variables  at  a  given  time  level  were  found  by  repeated 
iterations  of  the  governing  equations. 

Except  for  the  preset  isothermal  boundary  conditions,  the 
dependent  variables  were  first  initialized  by  setting  them  all  to  zero 
(However,  as  will  be  explained  later,  initial  conditions  other  than 
zero  were  used  when  searching  for  hysteresis  loops.) 

For  the  first  several  time  steps,  the  iterative  solution  process 
was  initiated  as  follows.  First,  the  E  correction  terms  were  set  to 
zero  and  a  converged  first-order  upwind  differenced  solution  was 
obtained.  This  converged  result  was  then  used  as  an  initial  condition 
for  obtaining  a  solution  to  the  corrected  central  differenced  scheme. 
This  two-step  method  allowed  for  an  enhanced  convergence  rate  when  the 
flow  field  was  diffusion-dominated--usual ly  associated  with  the  first 
few  time  steps  of  a  developing  flow.  However,  for  subsequent  time 
steps,  when  convective  effects  became  more  significant,  the  correction 
terms  were  always  retained  and  previously  converged  old-time  second- 


order  accurate  results  were  used  as  an  initial  start-up  for  the  next 
time  level.  The  advantages  of  this  procedure  are  more  fully  explained 
and  discussed  in  Kassinos  (1986).  Using  this  approach  allowed  for 
numerical  solutions  to  be  efficiently  computed  to  any  arbitrary  time. 

5.2.1.  Iteration  sequence  and  convergence  criterion 

The  governing  finite-difference  equations  (both  the  2-D  Navier- 
Stokes  and  the  boundary-layer  equations)  were  numerically  iterated  in 
the  following  sequence: 

i.  energy  equation 

ii.  vorticity  equation 

iii.  stream  function  equation. 

This  sequence  was  repeated  until  the  iterations  converged  to  within 
a  prescribed  tolerance.  The  maximum  modulus  of  the  difference  between 
temperature,  vorticity  and  stream  function  for  two  successive  iteration 
values  was  used  as  the  relative  convergence  constraint,  namely, 


max 


m+1 ,  n+1  m,  n+1 


(*m+1 •  n+1 ) 


max 


<  1  x  10 


-6 


(5.16) 


where  4>  represents  temperature,  vorticity  or  stream  function.  The 
iterative  process  was  terminated  and  the  numerical  solution  was 
considered  converged  when  the  above  criterion  was  satisfied.  In 
Eq.  (5.16),  the  m  referred  to  the  iteration  number,  while  the  time 
level  index  was  denoted  by  n.  The  vorticity  iterations  usually 
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3 

converged  the  slowest  (about  10  times  slower  when  compared  to 
temperature),  especially  for  highly  convective  or  multicellular-type 
flows. 


5.2.2.  Relaxation  parameters 


Two  independent  relaxation  parameters,  and  ^ >  were 
incorporated  into  the  equations  in  order  to  gain  some  control  on  the 
rate  of  convergence  of  the  iterations.  The  first  relaxation  parameter. 


,  was  used  only  with  the  vorticity  boundary  conditions  for  both  the 


2-D  Navier-Stokes  and  the  boundary-layer  equations.  Its  use  on  the 


inner  cylinder  boundary  can  be  illustrated  as  follows  (for  the 


Navier-Stokes  equations): 


wf+]’  n+1  =  {2f^ ’ jn+1/ (Gh1 )2} 

+  (1  -  ^)w^’jn+1  along  r  =  0  .  (5.17) 

If  <1,  then  the  m  +  1  iterated  value  of  w  becomes  a  weighted 
average  of  the  newly  computed  value  and  the  old  m  iterated  value  of  w. 
Generally,  .1  <_  ^  £  .5  was  used  to  help  stabilize  the  numerical 
computations.  When  multicellular  flow  was  encountered,  the  relaxation 
parameter  was  usually  reduced  to  the  smaller  values.  The  above 
procedure  was  also  used  for  the  boundary-layer  equations. 

Since  the  central  differenced  correction  terms  of  the  2-D 
Navier-Stokes  equations  were  based  on  the  old-time  level  (En),  an 
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underrelaxation  procedure  was  not  needed  to  aid  in  computational 
stability.  However,  for  the  finite-Prandtl  and  the  zero-Prandtl 
number  boundary-layer  equations,  a  second  relaxation  parameter,  , 
was  used  with  the  second-order  upwind  (or  central)  differenced 
correction  terms  which  were  based  on  the  new-time  level  (En+^).  The 
interior  point  equations  for  the  dependent  variables  <j>  =  T  and  w  for  the 
finite-Prandtl  number  equations,  and  4>  =  W  for  the  zero-Prandtl  number 
equations,  were  computed  according  to: 


m+1 ,  n+1 
<t>0 


,  r  m+1 ,  n+1  .  r  .m+1 ,  n+1  .  cm+l ,  n+1 ,  /r 
+  U<{>7  +  +  S  }/Cr 


+  {ft2  Em+1 ’  n+1  +  (1  -  ft2)  Em’  n+1}/CQ  .  (5.18) 


In  Eq.  (5.18),  the  Em  term  uses  4>q,  <f>.j ,  and  4>2  evaluated  at  the  m  -  1 
iteration  and  n  +  1  time  level,  together  with  <}>3  and  <p4  evaluated  at 
the  m  iteration  and  n  +  1  time  level.  Also,  <J>q  does  not  require  an 
iteration  index  since  it  represents  the  previously  converged  value 
at  the  old-time  level.  For  il2  <  1 ,  the  differencing  of  the  correction 
term  becomes  a  weighted  average  of  Gauss-Seidel  (current  level)  and 
Jacobi  (previous  level)  iterations.  Generally,  .1  <  ^  -  *5  was  used 
in  the  numerical  computations.  For  the  highly  convective  flows,  the 
relaxation  parameter  was  decreased  to  .2  or  .1. 


128 


5.2.3.  Multicellular  flow  determination 

In  the  2-D  Navier-Stokes  equations,  multicellular  flow  near  the 
top  of  the  annulus  was  investigated  for  Pr  =  .706.  When  searching 
for  the  transitional  Rayleigh  numbers  associated  with  particular  narrow 
gap  widths  (numbers  which  characterize  the  transition  from  bicellular 
to  multicellular  flow),  "zero"  initial  conditions  were  used  for  only 
the  first  pretransitional  Rayleigh  number  case.  Then,  for  each 
successive  slightly  higher  Rayleigh  number  (Ra)  run,  the  previous 
converged  Ra  solution  was  used  to  start  the  new  solution.  This 
procedure  was  repeated  until  at  some  Ra  a  transition  to  multicells 
occurred.  In  order  to  find  hysteresis  loops,  prescribed  multicellular 
initial  conditions  (usually  associated  with  the  multicellular  flow 
field  characterized  by  the  transitional  Rayleigh  number)  were  used  for 
gradually  backtracking  the  Ra  (also  in  a  successive  manner)  until  the 
bicellular  solution  again  became  apparent.  Different  solutions, 
dependent  upon  the  initial  conditions,  are  possible  as  a  result  of 
nonuniqueness  inherent  in  nonlinear-type  problems. 

For  the  boundary-layer  equations,  a  multicellular  flow  occurring 
near  the  vertical  portions  of  a  narrow  annulus  was  investigated  for 
the  small-Prandtl  number  cases.  The  values  of  G  or  G  (which  are 
proportional  to  Rayleigh  number)  corresponding  to  multicellular 
transition  were  determined  in  a  manner  similar  to  that  just  described, 

^  /V 

but  with  G  or  G  replacing  Ra  in  the  above  procedure. 


The  majority  of  the  calculations  were  carried  out  with  31  radial 
nodes  and  102  angular  nodes  for  either  the  0  -  it  or  the  0  -  2ir  test 
cases.  Using  the  variable-mesh  algorithm,  the  mesh  spacing  was 
slightly  reduced  (generally  by  a  factor  of  1.5  for  the  largest  to 
smallest  increment  size)  near  the  inner  and  outer  boundary-layers  and 
also  in  the  thermal  plume  region  (for  air)  near  the  top  of  the  inner 
cyl inder. 

For  the  2-D  Navier-Stokes  equations,  steady-state  pretransitional 
solutions  were  obtained  by  using  a  time  step  of  5  x  10"^.  Typically, 
300  time  steps  (t  =  .15)  were  required  to  reach  a  convincing  steady- 
state  result.  The  program  was  generally  more  stable  with  smaller  time 
steps.  This  was  probably  due  to  increased  diagonal  dominance,  since 
the  Cq  coefficient  gets  larger  in  magnitude  as  x  decreases  (see 
Eq.  5.6b).  For  larger  time  steps  (between  .001  and  .01),  the 
solution  would  often  diverge  as  the  Rayleigh  number  increased;  this  was 
especially  true  for  the  multicellular  flow  calculations.  Thus,  a 
small  constant  time  step  was  used  for  most  of  the  computations. 

A  typical  pretransitional  steady-state  result  took  approximately 

5-10  CPU  hours  on  the  Perkin-Elmer  3242  computer  system  (where  1  unit 

of  CRAY1  time  corresponds  to  approximately  50-60  units  of  Perkin-Elmer 

time),  or  an  average  of  1-2  minutes  per  time  step.  The  average  number 

of  iterations  per  time  step  was  approximately  20-30.  However,  when 

the  transition  to  multicells  occurred  (two-to-four  or  two-to-six 

-4 

cells),  the  time  step  was  maintained  at  5  x  10  ,  but  the  amount  of 


CPU  time  needed  to  achieve  steady-state  (t  =  .15)  increased  to 
approximately  30  hours  or  an  average  of  six  minutes  (for  approximately 
100  iterations)  per  time  step. 

For  two  of  the  narrow  gaps  (G  =  .100  and  G  =  .200)  associated  with 
hysteresis  behavior  (Pr  =  .706),  the  complete  annular  gap  was 
determined  numerically.  This  was  done  in  order  to  ascertain  if 
symmetric  cell  development  ensued  near  the  top  of  the  annulus,  as 
initially  expected.  In  both  of  these  cases,  symmetric  cellular 
patterns  resulted.  Based  on  this  information,  all  other  narrow  gaps 
for  air  were  numerically  studied  by  considering  only  the  half-annulus 
with  symmetry  boundary  conditions  imposed  at  'P  =  0  and  180°.  A 
31  x  102  mesh  was  again  used  for  these  symmetric  cases,  thus  creating 
an  extremely  fine  grid  for  the  half-plane  geometry.  Aside  from  G  =  .100 
and  .200,  the  transitional  Rayleigh  numbers  predicted  for  the  narrow 
gap  widths  were  calculated  by  assuming  symmetry  about  the  vertical 
center-line.  Note  that  for  G  =  .100  and  .200  (where  the  complete 
annulus  was  used),  a  31  x  102  mesh  was  still  employed,  but  the  grid 
points  were  significantly  concentrated  near  the  top  portion  of  the 
annular  gap  (by  factors  of  2.1  and  1.7,  respectively)  in  order  to  achieve 
the  necessary  accuracy  to  resolve  the  synmetric  multicellular  flow 
field.  For  the  G  =  .100  case,  an  erroneous  asymmetric  cellular  pattern 
actually  resulted  when  using  a  mesh  with  nodes  concentrated  near  the  top 
by  a  factor  of  1.7.  However,  a  six-cellular  symmetric  pattern  became 
evident  when  the  nodes  were  further  concentrated  to  a  factor  of  2.1. 


For  the  Pr  -*  0  boundary-layer  equations,  symmetry  about  the 
vertical  center-line  was  assumed  for  a  31  x  102  size  mesh  (this 


assumption  is  discussed  in  detail  in  the  latter  part  of  this  subsection). 
Prior  to  the  point  of  multicellular  instability  in  the  vertical  portion 
of  the  annulus  (G  <  5.1),  steady-state  solutions  were  easily  achieved 

A 

within  50  time  steps  for  a  At  of  1.0.  This  size  time  step  was  also 

~  1/2 

relatively  small  since  t  =  (Ra/Pr)  t.  Hence,  as  the  Rayleigh 

-1/2 

number  increased,  t,  in  turn,  decreased  in  proportion  to  Ra  .  The 
amount  of  CPU-time  required  to  reach  the  steady-state  condition  was 
approximately  three  hours  or  an  average  of  3.6  minutes  (for 
approximately  80  iterations)  per  time  step.  However,  when  the 

A 

multicellular  flow  instability  set  in  at  G  =  5.2,  the  steady-state 
condition  could  no  longer  be  achieved.  In  order  to  obtain  an  efficient 
and  time-accurate  response  of  the  unsteady  flow  field,  five  different  time 
steps  were  considered:  At  =  .1,  .25,  .5,  1.0  and  2.0.  The  first  four 
time  steps  provided  time-responses  of  the  stream  function  at  r  =  .5/^ 

=  90°  to  within  1  percent  error  of  each  other.  Also,  all  four  of  these 
time  steps  were  able  to  resolve  the  small  but  sudden  variations  of 
the  stream  function  with  time.  For  At  =  2.0,  a  notable  difference  in 
the  time-response  was  observed.  That  is,  the  small  sudden  variations 
of  the  stream  function  with  time  could  no  longer  be  resolved.  Values 
of  the  stream  function  were  in  error  by  approximately  4-5  percent  when 
compared  to  those  corresponding  to  At  =  .10.  These  differences  are 
most  likely  due  to  increased  time  truncation  error  when  using  At  =  2.0. 

In  addition,  for  At  =  2.0,  the  number  of  iterations  required  for 
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convergence  at  a  particular  time  level  was  disproportionately  larger 

A  yv 

when  compared  to  that  for  At  =  1.0.  For  example,  an  increase  in  t 
from  35.0  to  36.0  (At  =  1.0)  took  approximately  550  iterations, 

A  /V 

whereas  an  increase  in  t  from  35.0  to  37.0  (At  =  2.0)  required  about 
1375  iterations  for  convergence.  This  was  probably  due  to  a  slight 
loss  in  numerical  stability  as  a  result  of  a  decrease  in  the  CQ 

/V 

coefficient  (see  Eq.  5.10b)  when  At  increases  from  1.0  to  2.0. 

Thus,  due  to  computational  time  and  accuracy  considerations,  the 
time  step  of  1.0  was  chosen  to  numerically  predict  the  time-response 

/V  A 

of  the  multicellular  flow  at  G  =  5.2.  Therefore,  for  values  of  G 
slightly  greater  than  5.2,  a  time  step  of  1.0  or  less  should  be  employed 
For  At  =  1.0,  the  average  CPU-time  per  time  step  amounted  to  12-14 

A 

minutes  (for  approximately  300  iterations)  with  G  =  5.2  and  a  half- 
annular  mesh  of  31  x  102  nodes.  The  number  of  nodes  in  the  angular 
direction  was  based  on  the  maximum  possible  to  properly  capture  cellular 
development  while  considering  computer  power  and  time  limitations.  In 
Section  6.2,  a  grid  study  for  four  different  size  meshes  showed  that 
the  31  x  102  node  mesh  was  fine  enough  to  resolve  the  same  number  of 
initial  cells  (seven)  as  the  finer  mesh  of  41  x  132.  In  further 
support  of  using  the  31  x  102  mesh,  Lee  and  Korpela  (1983)  indicated 
that  10  grid  points  per  cell  were  adequate  for  properly  resolving 
cellular  development  in  multicellular  natural  convective  flow  between 
long  and  narrow  vertical  slots.  Therefore,  since  the  instability  at 
G  =  5.2  resulted  in  7-8  cells  per  half-annulus,  102  nodes  in  the 
angular  direction  satisfied  their  particular  recommendation  for 
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avoiding  significant  resolution  problems. 

The  use  of  symmetry  boundary  conditions  with  regard  to  unsteady 
cellular  behavior  in  the  half-annulus  can  be  supported  by  the  nature 
of  the  cellular  activity  that  occurred.  The  strongest  cellular 
behavior  took  place  near  the  90°-point  vertical  section  of  the  annulus, 
whereas  near  ^  =  0°  and  ^  =  180°,  the  flow  field  remained  quiet  and 
undisturbed,  preserving  the  end  portions  of  the  familiar  steady-state 
kidney-shaped  flow  pattern.  In  support  of  this  claim,  the  analytical 
steady-state  result  derived  in  Chapter  4  (for  the  Pr  -*■  0  limit) 
demonstrated  that  the  steady-state  solutions  were  completely  symmetric 
about  '!'  =  0  and  180°.  Moreover,  the  cellular  development  in  the 
vertical  section  of  the  half-annulus  was  not  counter-rotating,  hence 
the  fluid  was  retained  within  the  outer  kidney-shaped  streamline.  This 
tends  to  indicate  that  the  cellular  flow  in  each  half-annulus  develops 
simultaneously  and  independently  of  each  other.  Also,  when  negative 
'P  is  substituted  into  the  Pr  -*■  0  governing  equations,  the  original 
form  of  the  equations  is  maintained--again  implying  symmetry  about 
=  0  and  180°. 

However,  in  order  to  properly  verify  this  salient  symmetry 
assumption,  the  symmetry  conditions  were  relaxed  and  the  flow  in  the 
complete  annular  gap  was  determined.  For  the  complete  annulus  (0  -  2t), 
a  31  x  202  mesh  was  used  and  the  initial  time-response  behavior  of 

A 

the  unsteady  cellular  motion  was  obtained  with  At  =  1.0.  The  initial 
conditions  used  to  start  the  run  were  provided  by  the  half-annular 
7-cell  solution  at  t  =  35.0  (a  mirror-image  of  this  solution  was  used 


to  incorporate  the  whole  annulus).  The  initial  time-response  behavior 
(the  first  832  time  steps  corresponding  to  one  cycle)  of  the  0  -  2tt 
run  reproduced  that  of  the  0  -  tt  run  within  one  percent  error  (see 
Figure  6.14).  The  important  changes  of  the  cellular  pattern  with 
time  were  also  duplicated.  For  the  complete  annulus,  the  average 
CPU-time  per  time  step  was  about  26  minutes  (for  approximately  400 
iterations). 

Based  on  this  single  test  case  and  preceding  relevant  supportive 
discussion,  the  half-annular  symmetry  assumption  was  used  throughout 
this  research  to  predict  the  unsteady  cellular  behavior  associated  with 
the  vertical  portions  of  a  narrow  annulus  (for  Pr  -*■  0).  Proceeding 
with  this  symmetry  assumption  substantially  reduced  the  CPU-time 
necessary  for  a  converged  solution,  thus  rendering  the  problem  at 
hand  feasible  to  study  on  the  Perkin-Elmer  3242  computer  system. 

The  numerical  handling  of  the  finite-Prandtl  number  boundary-layer 
equations  was  very  similar  to  that  previously  described  for  the  Pr  -*■  0 
equations.  These  equations  were  also  solved  assuming  symmetry  using 
a  hal f-aunular  mesh  of  31  x  102  nodes  with  At  =  1.0.  The  finite- 
Prandtl  number  equations  were  usually  solved  numerically  to  determine 
the  transitional  values  of  G,  for  various  low-Prandtl  number  fluids, 
corresponding  to  the  onset  of  multicellular  instability.  These 
equations  were  not  used  to  study  any  type  of  unsteady  behavior  relating 
to  a  finite-Prandtl  number  hydrodynamic  instability. 

Upon  solving  the  equations,  it  was  noticed  that  as  Prandtl  number 
increased,  so  did  the  computation  time  necessary  to  reach  convergence 


at  a  particular  time  level.  This  was  because  the  governing  equations 
were  nondimensional ized  with  respect  to  momentum  diffusivity  as  opposed 
to  thermal  diffusivity,  thus  allowing  the  energy  equation  to  become 
diffusion  dominated  for  small  Prandtl  numbers.  To  reverse  this 
situation,  the  governing  equations  must  be  recast  to  allow  the  vorticity 
equation  to  become  diffusion  dominated  as  Prandtl  number  increases. 

Lastly,  it  should  be  noted  that  for  both  the  2-D  Navier-Stokes 
and  finite-Prandtl  number  boundary-layer  equations,  another  important 
measure  of  numerical  accuracy  can  be  determined  for  the  steady-state 
cases.  That  is,  the  average  inner  Nusselt  number  (Nu^)  should  equal 
the  average  outer  Nusselt  number  (Nuq)  at  steady-state.  For  all  of  the 
steady-state  test  results,  the  equality  was  in  error  by  less  than 
.1  percent.  For  the  pretransitional  results  of  the  2-D  Navier-Stokes 
equations,  Nuq  was  slightly  less  than  Nu^ .  However,  for  the 
multicellular  flows,  Nun  was  slightly  greater  than  Nu.. 


6.  RESULTS  AND  DISCUSSION 

This  chapter  consists  of  two  sections.  The  first  essentially 
describes  the  hysteresis  behavior  associated  with  steady  multicellular 
flow  occurring  in  narrow  gaps  for  air.  This  type  of  multicellular  flow 
apparently  originates  because  of  thermal  instability,  and  develops  near 
the  top  portion  of  a  narrow  horizontal  annulus. 

The  second  section  mainly  describes  an  unsteady  multicellular  flow 
that  develops  near  the  vertical  portions  of  a  narrow  horizontal  annulus 
for  smal 1 -Prandtl  number  fluids.  This  particular  type  of  instability 
appears  to  be  hydrodynamic  in  origin. 

Wherever  possible,  results  pertinent  to  this  study  will  be  compared 
to  related  analytical,  numerical  and/or  experimental  work. 

6.1.  Thermal  Instability 
6.1.1.  Hysteresis  behavior 

As  described  in  Chapter  5,  the  complete  annular  flow  field  (0-2tt) 
was  determined  numerically  when  searching  for  hysteresis  behavior.  A 
31  x  102  mesh  was  employed  with  the  grid  points  concentrated  near  the 
top  horizontal  portion  of  the  annulus  -  where  the  multicellular  flow 
field  tended  to  develop  for  air. 

The  average  Nusselt  number  and  shear-stress  results  are  presented 


in  Figures  6.1a,  6.1b  and  6.2,  respectively,  for  a  gap  number  G  =  .200 
and  Pr  =  .706  (note  that  these  average  dimensionless  parameters  are 
defined  in  Appendix  D).  The  transition  of  the  flow  field  from  two 


Hysteresis  behavior  of  the  mean  Nusselt  number  variation  with  Rayleigh  number 


Hysteresis  loop  of  Figure  6.1a  with  an  extended 
number  variation 


Hysteresis  behavior  of  the  average 
Rayleigh  number  for  G  =  .200  and  P 


to  four  cells  occurred  abruptly  at  a  Rayleigh  number  of  approximately 
351,000  (or  Ra^_a  =  2808).  Associated  with  this  transition  was  a 
sudden  rise  ir.  the  average  Nusselt  number  (see  Figure  6.1a)  due  mainly 
to  mcr?  efficient  fluid  mixing  created  by  the  counter-rotating 
secondary  cellular  motion.  This  is  analogous  to  the  increase  in  heat 
transfer  experienced  when  a  flow  makes  the  transition  from  laminar  to 
turbulent  fluid  motion.  Note  that  the  above  four-cell  transition  does 
not  indicate  the  onset  of  turbulence,  but  rather  signifies  a  more  complex 
laminar  multicellular  flow  condition.  Also,  because  of  the  two  smaller 
counter-rotating  cells  initiated  near  the  top  of  the  annulus,  the  local 
shear-stress  in  this  region  reversed  its  sign.  Hence,  a  sudden 
decrease  in  the  average  inner  shear-stress  (for  the  half-annulus) 
occurred  at  the  transition  from  two  to  four  cells  (see  Figure  6.2). 

This  particular  phenomena  bears  resemblance  to  a  similar  occurrence 
in  the  Taylor-Instability  problem  of  Coles  (1965),  for  fluid  flow  between 
vertical  concentric  narrow  rotating  cylinders.  In  this  problem,  a 
hysteresis  effect  of  the  torque  on  the  inner  cylinder  (which  is 
directly  proportional  to  shear-stress)  was  observed  by  Eagles  (1974). 

He  showed  that  at  a  certain  Taylor  number,  the  Taylor-vortex  flow 
would  change  to  a  wavy-vortex  flow,  and  associated  with  this  change 
was  a  sudden  lowering  of  the  inner  cylinder  torque  (as  also  discussed 
in  Stuart,  1960),  although  the  Taylor  number  had  actually  increased 
in  magnitude. 

Both  the  Nusselt  number  and  the  shear-stress  data  exhibited  a 
hysteresis  behavior  as  the  Rayleigh  number  was  gradually  increased 


and  then  decreased  past  the  transition  point.  For  Ra  >  350,000,  only 
a  four-cell  solution  prevailed,  while  for  Ra  <  285,000,  only  a  two-cell 
solution  prevailed.  However,  between  these  two  limits,  both  types  of 
solutions  were  possible,  thus  forming  a  hysteresis  loop.  The  numerical 
approach  used  to  determine  the  transitional  Rayleigh  numbers  and  related 
hysteresis  behavior  was  described  in  detail  in  Chapter  5. 

Powe  et  al_.  (1971)  numerically  estimated  a  transitional  Rayleigh 
number  of  452,000  for  G  =  .200  based  upon  the  point  at  which  the  stream 
function  changed  its  sign.  However,  they  could  not  fully  resolve  the 
smaller  counter-rotating  cell  with  their  particular  method  (see 
Chapter  2).  They  used  a  uniform  increment  mesh  with  35  angular  and 
15  radial  nodes  for  the  half-annulus.  The  fact  that  they  used  a 
relatively  coarse  mesh  in  the  angular  direction  may  have  caused 
resolution  problems  which  delayed  their  transitional  Rayleigh  number 
to  a  value  approximately  30  percent  higher  than  the  351,000  value 
predicted  in  this  study.  A  similar  four-cellular  flow  pattern  was 
experimentally  observed  by  Bishop  e_t  al_.  (1964)  in  a  concentric  spherical 
geometry.  For  G  =  .19,  they  obtained  a  transitional  Rayleigh  number 
(based  on  gap  width)  of  about  3600,  which  corresponds  to  approximately 
Ra  =  525,000. 

These  earlier  results  agree  with  the  present  work  in  that  a 
four-cell  pattern  predominates  immediately  after  the  thermal  instability 
sets  in  for  G  =  .200. 

To  improve  the  numerical  prediction  of  the  corresponding 
transitional  Rayleigh  number,  grid  nodes  should  be  concentrated  in  the 


thermal  plume  region  where  the  smaller  counter-rotating  cells  develop, 
as  was  properly  done  in  this  study  (see  Section  5.2.4  and  Appendix  E 
for  further  discussion). 

Streamline  plots  of  the  cellular  flow  field  development  for 
G  =  .200  are  depicted  in  Figures  6.3a-6.3c.  Figure  6.3a  corresponds 
to  a  pretransitional  result,  Ra  =  350,000,  and  demonstrates  the  usual 
kidney-shaped  flow  pattern.  Figure  6.3b  shows  the  initial  development 
of  the  two  smaller  counter-rotating  cells  for  Ra  =  351,000.  Figure  6.3c 
depicts  the  secondary  flow  in  a  much  stronger  state  beyond  transition, 
at  Ra  =  900,000.  A  typical  multicellular  temperature  profile  is  also 
given  in  Figure  6.3d  for  Ra  =  900,000.  The  shape  of  the  isotherms  was 
responsive  to  the  fluid  motion  for  Pr  =  .706,  resulting  in  an  inverted 
thermal  plume  near  the  top  portion  of  the  annulus.  In  Figures  6.3a-d, 
the  diameter  ratio  is  shown  as  1.5  times  its  actual  size  for  ease  of 
viewing. 

Similarly,  for  G  =  .100,  an  abrupt  rise  in  heat  transfer  occurred 
as  the  flow  field  made  the  transition  from  two  to  six  cells  at 
RaTR  =  2,841,000  (or  Ra^_a  =  2841),  as  shown  in  Figure  6.4.  The 
hysteresis  loop  for  this  case  expanded  from  Ra  =  2,570,000  to 
Ra  =  2,840,000,  a  change  of  approximately  270,000.  Below  this  range, 
only  the  two-cell  solution  prevailed,  while  above  this  range,  only  the 
six-cell  solution  prevailed.  For  G  =  .200,  the  Rayleigh  number  change 
in  the  hysteresis  loop  amounted  to  approximately  65,000,  considerably 
smaller  than  for  G  =  .100.  Furthermore,  the  rise  in  the  average  Nusselt 
number  associated  with  the  G  =  .200  multicellular  transition  was  greater 


Figure  6.3b.  Streamlines  of  an  initial  four-cellular  flow  field 
at  Ra  =  351,000  for  G  =  .200  and  Pr  =  .706 
(diameter  ratio  shown  as  1.5  times  actual  size) 


3c.  Streamlines  of  a  more  developed  four-cellular  flow 
field  at  Ra  =  900,000  for  G  =  .200  and  Pr  =  .706 
(diameter  ratio  shown  as  1.5  times  actual  size) 


gure  6.3d.  Isotherms  of  a  more  developed  four-cellular  flow 
field  at  Ra  =  900,000  for  G  =  .200  and  Pr  =  .706 
(diameter  ratio  shown  as  1.5  times  actual  size) 


Hysteresis  behavior  of  the  mean  Nusselt  number  variation  with  Rayleigh  number 


than  for  G  =  .100  by  an  amount  of  approximately  .052,  or  about  a  50 
percent  increase. 

The  above  occurrences  may  be  explained  as  follows.  As  the  gap 
decreases  in  size,  both  the  viscous  effects  and  the  restraining 
influence  of  the  geometry  become  more  significant  for  a  given  value  of 
Rayleigh  number,  thus  suppressing  natural  convective  fluid  motion  within 
the  annulus.  This  has  several  ramifications.  First,  because  of  the 
fluid  motion  hindrance,  a  higher  transitional  Rayleigh  number  is  needed 
to  trigger  the  instability  for  G  =  .100,  compared  to  that  for  G  =  .200. 
Secondly,  once  the  multicellular  instability  results,  related  hysteresis 
effects  persist  longer  for  G  =  .100  due  to  suppressed  fluid  motion 
changes.  Hence,  a  greater  Rayleigh  number  range  for  the  hysteresis 
loop  of  this  gap  occurs.  Thirdly,  suppressed  convective  effects  for 
G  =  .100  may  also  cause  a  less  pronounced  rise  in  the  average  Nusselt 
number  associated  with  multicellular  transition,  compared  to  that  for 
G  =  .200.  Lastly,  as  the  gap  decreases,  curvature  effects  diminish, 
resulting  in  a  better  comparison  to  the  Benard  type  of  instability 
between  horizontal  parallel  flat  plates.  This,  coupled  with  increased 
viscous  effects,  may  cause  the  number  of  cells  occurring  at  transition 
to  increase,  i.e.,  a  two  to  six  cell  transition  for  G  =  .100  as  opposed 
to  the  two  to  four  cell  transition  with  G  =  .200. 

The  results  of  this  study  seem  to  indicate  that  for  G  <  .175  +  .015 
an  initial  transition  from  two  to  six  cells  occurs,  rather  than  the 
two  to  four  cell  transition  for  G  =  .200.  This  statement  is  supported 
by  the  numerical  results  of  Rao  e_t  al_.  (1985),  where  they  graphically 


depicted  a  six-cell  flow  pattern  at  a  Rayleigh  number  of  approximately 
750,000  for  6  =  .175.  However,  they  did  not  give  the  Rayleigh  number 
corresponding  to  the  point  of  transition  for  this  particular  gap  size. 

In  this  study,  for  G  =  .175,  a  transition  to  six-cells  took  place  at 
Ra  =  450,000.  A  typical  six-cell  flow  pattern  with  related  isotherms, 
for  G  =  .100  and  Ra  =  3,600,000,  is  depicted  in  Figures  6.5  and  6.6, 
respectively.  Notice  the  upright  thermal  plume  (Figure  6.6)  associated 
with  the  six-cell  flow  field. 

For  G  =  .100,  a  hysteresis  behavior  was  also  present  for  the 
average  inner  shear-stress  versus  Rayleigh  number.  However,  these 
results  were  not  plotted  since  the  related  hysteresis  loop  was  very 
slight  and  almost  indistinguishable  in  a  graphical  sense--where  only 
a  .2  percent  increase  in  the  average  inner  shear-stress  occurred  at  the 
point  of  multicellular  transition  for  G  =  .100. 

Figures  6.7  and  6.8  indicate  the  effect  of  multicellular  flow  on 
the  angular  variation  of  the  inner  Nusselt  number  and  centerline  stream 
function,  respectively,  with  regard  to  G  =  .100  and  G  =  .200.  To 
understand  these  effects,  cell  rotation  must  be  taken  into  consideration. 
For  the  half-annulus  cellular  pattern  of  G  =  .200,  the  smaller  secondary 
cell  rotates  counter-clockwise  and  extends  between  165°  and  180°  for 
Ra  =  900,000  (see  Figure  6.3c).  This  rotation  is  responsible  for  the 
development  of  an  inverted  thermal  plume  near  the  top  portion  of  the 
annulus.  For  the  cellular  pattern  of  G  =  .100  at  Ra  =  3,600,000,  two 
secondary  cells  are  present  in  the  half-annulus.  One  occurs  between 
170°  and  180°  and  rotates  clockwise,  while  the  other  occurs  between 


Figure  6.5.  Streamlines  of  a  six-cellular  flow  field  at 
Ra  =  3,600,000  for  G  =  .100  and  Pr  =  .706 
(diameter  ratio  shown  as  1.5  times  actual  s 


Figure  6.7.  Angular  variation  of  the  inner  (r=0)  Nus 
.200  at  Ra  =  3,600,000  and  900,000,  respi 


165°  and  170°  and  rotates  counter-clockwise  (see  Figure  6.5). 


Due  to  the  inverted  thermal  plume  for  G  =  .200,  Figure  6.7  shows 
a  local  maximum  in  heat  transfer  on  the  inner  wall  at  4>  =  180°  and  a 
local  minimum  near  4*  =  165°.  In  contrast,  the  six-cell  flow  pattern 
of  G  =  .100  shows  a  local  maximum  in  heat  transfer  near  4^  =  170°  due 
to  the  counter-clockwise  rotating  secondary  cell,  while  a  local  minimum 
results  at  4^  =  180°.  This  minimum  occurs  because  the  clockwise 
rotating  secondary  cell  causes  an  upright  thermal  plume  near  the  top 
of  the  inner  cylinder.  Also,  smaller  local  maximum  and  minimum  in  heat 
transfer  result  near  4^  =  155°  and  165°,  respectively,  due  to  the  thermal 
response  of  the  clockwise  rotating  primary  cell  (which  bottle-necks) 
near  4>  =  155°  coupled  with  the  counter-clockwise  rotating  secondary  cell. 

In  a  similar  manner,  the  two  types  of  multicellular  flow  also  cause 
a  departure  of  the  unicellular  sinusoidal -1 i ke  stream  function 
distribution  between  41  =  150°  and  180°.  This  departure  is  illustrated 
in  Figure  6.8  for  the  four  and  six-cell  flow  patterns. 

6.1.2.  Analytical  comparison 

The  analytical  perturbation  expressions  derived  in  Chapter  4, 
from  the  fini te-Prandtl  number  boundary- layer  equations,  were  employed 
as  convenient  checks  to  the  pretransitional  results  generated  in  the 
preceding  numerical  analysis. 

The  forthcoming  comparison  has  a  dual  purpose:  first,  it  is 
valuable  in  the  sense  that  the  amount  of  available  data  pertaining  to 
local  flow  changes  in  the  narrow-type  gaps  is  extremely  scarce  for  the 
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high  Rayleigh  number  range;  and  second,  a  valid  comparison  will  serve 
to  support  the  numerical  method  used  in  this  study,  which  in  turn, 
supports  the  asymptotic  analysis  described  in  Chapter  4. 

The  three-term  inner  Nusselt  number  perturbation  expression  given 
in  Eq.  (4.54),  together  with  the  three-term  expressions  for  vorticity 
and  stream  function  given  in  Eq.  (4.29),  will  be  used  to  generate 
comparative  data  for  the  pretransitional  2-D  Navier-Stokes  numerical 
results.  Comparisons  for  the  inner  (r  =  0)  Nusselt  number,  the  inner- 
wall  (r  =  0)  vorticity  and  the  center-line  (r  =  .5)  stream  function  are 
provided  in  Figures  6.9-6.11,  respectively. 

Figure  6.9  depicts  the  angular  variation  of  the  inner  Nusselt 
number  using  the  analytical  expression  (4.54)  with  G  =  4.0  and  Pr  =  .706. 
It  is  worth  noting  that  G  =  4.0  does  not  correspond  to  a  conduction- 
dominated  flow  (as  with  G  =  1.0  or  2.0),  but  rather  a  flow  subject  to 
significant  convective  effects.  Recalling  that  G  =  Ra^G,  four 
different  combinations  of  G  and  Ra  corresponding  to  G  =  4.0  will  be 
used  as  numerical  comparisons  to  the  analytical  result.  These 
combinations  are  G  =  .200  (Ra  =  160,000),  G  =  .150  (Ra  =  505,680), 

G  =  .100  (Ra  =  2,560,000)  and  G  =  .025  (Ra  =  6.554  x  108).  Related 
results  are  also  displayed  in  Figure  6.9. 

The  best  comparison  to  the  analytical  result  corresponded  to 
G  =  .025.  This  was  expected  since  the  asymptotic  analysis  was  tailored 
for  the  smal 1 -gap/high  Rayleigh  number  flow  regime.  Therefore,  the  2-D 
Navier-Stokes  solutions  do  indeed  converge  toward  the  perturbation 
solutions  as  G  -*■  0.  In  the  comparisons  that  follow,  relative  errors 


of  the  2-D  Navier-Stokes  results  with  respect  to  the  analytical  results 
are  calculated.  For  G  =  .025,  a  maximum  error  of  1.5  percent  resulted 
near  the  top  of  the  inner  cylinder,  ^  =  180°.  The  three  larger  size 
gaps  also  compared  favorably.  For  G  =  .200,  .150  and  .100,  maximum 
respective  errors  of  11.5,  8.5  and  6.0  percent  occurred,  also  at  'P  =  180° 

In  Figure  6.10,  a  comparison  similar  to  Figure  6.9  is  presented. 

This  new  comparison  regards  the  angular  variation  of  the  inner-wall 
vorticity  for  G  =  4.0  and  Pr  =  .706.  For  G  =  .200,  .100  and  .025, 
maximum  respective  errors  of  7.1,  3.5  and  1.0  percent  resulted  near 
^  =  90°,  or  approximately  the  point  of  maximum  fluid  velocity. 

Lastly,  Figure  6.11  depicts  the  angular  variation  of  the  center-line 
stream  function  for  G  =  4.0  and  Pr  =  .706.  This  particular  comparison 
was  in  very  close  agreement  to  the  2-D  Navier-Stokes  results,  with 
maximum  respective  errors  of  2.0,  1.5  and  1.0  percent  for  G  =  .200, 

.100  and  .025. 

Although  the  inner  Nusselt  number  (which  is  proportional  to  the 
radial  temperature  gradient)  is  generally  considered  the  least  sensitive 
parameter  of  the  three,  it  did  not  agree  as  closely  as  did  the 
inner-wall  vorticity  and  the  center-line  stream  function  comparisons. 

Most  probably,  this  agreement  would  improve  if  more  expansion  terms 
are  incorporated  into  the  analytical  expression  for  the  inner  Nusselt 
number.  But  overall,  all  three  analytical  comparisons  agreed  rather 
well  with  the  2-D  Navier-Stokes  results,  hence  supporting  both  the 
numerical  and  analytical  work  set  forth  in  this  thesis. 


gure  6.10.  Angular  variation  of  the  inner-wall  (r=0)  vorticity  (7  Ra  '  )  for 

G  =  G  Ra1//4  =  4.0  and  Pr  =  .706:  2-D  Navier-Stokes/analytical  comparison 
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In  summary,  the  first  three  terms  of  the  perturbation  solution 
appear  to  be  valid  to  within  approximately  10  percent  for  the 
pretransitional  cases  of  G  £  .200.  Thus,  quick  and  fairly  accurate 
approximate  solutions  to  the  2-D  Navier-Stokes  equations,  in  this 
particular  range  of  narrow-gap  widths,  can  be  obtained  from  the 
analytical  expressions  derived  in  this  study. 

6.1.3.  Small -gap  number  stability  curves 

As  discussed  in  Chapter  4,  previous  investigators  (Walton,  1980; 
Powe  et  aK  ,  1971;  Liu  et^  aj_. ,  1962)  speculated  that  as  the  narrow 
annular  gap  spacing  approaches  zero,  a  true  Benard-type  instability 
would  erupt,  corresponding  to  a  critical  Rayleigh  number  (Ra.  )  of 
approximately  1700.  Instead,  it  appears  from  Figure  6.12  that  the 
minimum  transitional  Rayleigh  number  approaches  2400  for  G  =  .175. 
Beyond  this  point,  the  transitional  Rayleigh  number  begins  to  increase 
substantially  with  decreasing  gap  number;  where  at  G  =  .050,  the 
thermal  instability  sets  in  at  approximately  Ra^  =  4560. 

Directly  related  to  the  stability  curve  of  Figure  6.12  is 
Figure  6.13,  which  was  obtained  from  the  relation  G  =  Ra^4G.  This 
new  curve  appears  monotonic  and  seems  to  indicate  that  for  very  small 
gaps,  the  instability  sets  in  at  a  near-constant  value  of  G  equal  to 
approximately  3.90. 

All  data  points  on  the  stability  curve  of  Figure  6.12  were 
calculated  assuming  vertical  symmetry  for  a  half-annular  mesh  of 
31  x  102  nodes,  except  for  G  ~  .100  and  .200  where  the  complete  annulus 


(0  -  2tt)  was  considered  when  investigating  possible  hysteresis  behavior 
The  vertical  symmetry  assumption  proved  to  be  valid,  as  was  discussed 
in  Section  5.2.4.  In  order  to  obtain  some  type  of  tolerance  criteria 
for  the  numerically  predicted  transitional  Rayleigh  numbers,  a  mesh 
resolution  study  was  performed  (see  Appendix  E)  for  three  different 
size  gap  widths.  Based  on  these  results,  it  appears  that  the  predicted 
transitional  Rayleigh  numbers  are  within  5  to  15  percent  uncertainty. 
This  error  range  seems  to  indicate  fairly  reasonable  accuracy  when 
compared  to  some  of  the  larger  discrepancies  reported  in  the  literature 
for  predicting  the  onset  of  secondary  motion  (Elder,  1965;  Powe  et  al . , 
1971;  Rao  et  ,  1985). 

Also  in  Appendix  E  is  a  comparison  of  the  cellular  development 
for  G  =  .200  when  using  either  a  first-order  upwind  or  a  second-order 
central-differencing  method.  This  comparison  demonstrated  that  even 
at  very  high  Rayleigh  numbers,  the  first-order  upwind  differencing 
scheme  could  not  capture  multicellular  flow  transition.  In  addition, 
this  appendix  compares  mean  Nusselt  numbers  predicted  numerically  for 
G  =  .200  to  that  estimated  by  the  correlation  of  Raithby  and 
Hollands  (1975),  which  was  valid  for  narrow  gaps  in  the  pretransi tional 
convecti ve-dominated  flow  regimes,  and  also  to  the  correlation  of 
Kuehn  and  Goldstein  (1980b). 

6.2.  Hydrodynamic  Instability 

Under  certain  conditions,  an  unsteady  multicellular  flow 
instability  originates  in  the  vertical  sections  of  two-dimensional 


narrow  horizontal  cylindrical  annuli.  This  particular  type  of 
multicellular  flow  was  captured  numerically  from  the  simplified 
governing  equations  derived  in  Section  4.5,  corresponding  to  the  high 
Rayleigh  number/narrow-gap/small -Prandtl  number  limiting  conditions. 

The  multicells  apparently  result  because  of  hydrodynamic 
instability,  and  are  analogous  to  those  investigated  by  other  researchers 
with  regard  to  the  vertical  slot  geometry,  as  discussed  in  detail  in 
Section  2.6. 

In  contrast  to  the  vertical  slot  studies,  the  multicellular  flow 
initiated  in  the  vertical  portions  of  a  horizontal  annulus  is  unsteady. 
Results  from  this  study  seem  to  indicate  that  the  initial  instability 
behaves  in  a  rather  periodic  manner.  That  is,  the  strength  of  the 
cells  appears  to  increase  and  decrease  with  time  in  a  cyclic  fashion. 
Closely  related  to  this  observation  is  the  phenomena  of  scrange- 
attractors.  In  some  systems  that  experience  a  multicellular-like  flow 
instability,  an  ordered  route  to  chaos  has  been  reported  (see  Section 
2.5  for  specific  references  and  related  discussion).  Briefly  though, 
a  typical  sequence  of  events  consists  of  the  system  first  behaving  in 
a  time-periodic  manner,  followed  by  a  cycle  of  periodic-doubling  until 
finally,  the  system  is  stressed  into  chaotic-like  nonperiodic  motion. 

This  present  flow  under  study  seems  to  have  the  potential  to 
follow  a  similar  ordered  path  to  chaos. 

The  following  subsections  deal  with  a  thorough  discussion  of  both 
the  numerical  and  analytical  results  relating  to  this  unique  type  of 
secondary-flow  instability. 


Assuming  symmetry  about  the  vertical  center-line  greatly  reduced 
the  amount  of  CPU-time  associated  with  the  numerical  computations. 

This  assumption  was  based  on  the  detailed  discussion  given  in  Section 
5.2.  The  above  symmetry  condition  was  proven  valid  in  Figure  6.14a. 

A 

Here,  the  first  832  time  steps  (for  At  =  1.0)  of  the  numerical  solution 
to  the  complete  annulus  (31  x  202  nodes)  produced  the  cellular  cycle  of 
the  half-annulus  (31  x  102  nodes)  almost  identically,  and  graphically, 
no  notable  differences  could  be  discerned.  Also,  for  both  the  full  and 

A 

half-annulus,  the  inner-wall  vorticity  at  t  -  832  were  compared  in 
Figure  6.14b.  The  two  cases  were  again  in  close  agreement,  to  within 
approximately  one  percent  error.  Note  that  these  curves  are  presented 
at  this  time  to  support  the  assumption  of  vertical  symmetry.  The  exact 
nature  and  meaning  of  the  cellular  behavior  described  by  these  curves 
will  be  explained  later  in  this  section.  Therefore,  for  the  results 
that  follow,  a  half-annulus  mesh  of  31  radial  nodes  and  102  angular 
nodes  was  employed.  In  general,  the  justification  for  this  size  mesh 
was  due  to  proper  resolution  of  the  cellular  structure  and  also  to 
computer  power  and  time  limitations.  These  considerations  are  again 
discussed  in  more  detail  in  Section  5.2. 

Steady-state  results  were  obtained  up  to  a  G  value  of  5.1.  Then, 
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at  G  =  5.2,  a  seven-cell  unsteady  instability  set  in  at  a  t  value  of 
approximately  30.0.  Note  that  the  converged  steady-state  results  for 
G  =  5.1  were  used  as  initial  conditions  to  start  the  G  =  5.2  case. 

Prior  to  the  seven-cell  development,  the  instability  first  occurs  near 
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ip  =  90°,  where  the  maximum  buoyancy  force  causes  the  most  rapid  rising 
and  falling  of  the  fluid,  as  illustrated  in  Figure  6.15a.  Also  in 
Figure  6.15b,  the  cellular  pattern  is  displayed  at  a  time  slightly 
beyond  the  initial  development  stage.  The  cellular  instability  then 
grew  until  the  relatively  stable  seven-cell  pattern  resulted  at 
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t  =  30.0.  The  initial  seven-cell  solution  encompassed  an  arc  of  about 
80  degrees  near  the  vertical  portion  of  the  annulus  (see  Figure  6.16). 

It  is  important  to  note  that  this  multicellular  flow  instability 
was  preserved  with  both  finer  and  coarser  mesh  sizes,  as  evidenced  by 
the  data  in  Table  6.1.  The  coarser  mesh,  31  x  72,  could  not  resolve  the 
smaller  cells  near  the  top  and  bottom  portion  of  the  vertical  section, 
thus  resulting  in  only  five  cells  for  the  initial  instability.  The 
other  three  mesh  sizes  indicate  that  the  instability  sets  in  as  a  six 
or  seven-cellular  flow  pattern.  This  result  was  in  close  agreement  with 
that  reported  by  Lee  and  Korpela  (1983)  for  their  vertical  slot  geometry 

Table  6.1.  Mesh  resolution  comparison 


Mesh 

size 

G 

transitional 

Geometry 

Cells 

Convergence 

constraint 

Precision 

31  x 

72 

5.7 

hal f-annul us 

5 

1  X 

IQ'6 

single 

31  x 

102 

5.2 

hal f-annul us 

7 

1  X 

10'6 

single 

41  x 

132 

4.9 

half-annulus 

7 

1  X 

10'8 

double 

31  X 

162 

5.2 

full-annulus 

6 

1  X 

10-6 

single 

172 


They  claimed  that  a  six-cellular  flow  pattern  initially  appeared  at 
Gr.  =  8000  for  a  Prandtl  number  very  near  zero  (see  Section  2.6). 

This  Grashof  number  of  8000  compared  favorably  with  the  transitional 
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value  of  this  study,  where  G  =  5.2  translated  to  a  Grashof  number  (based 
on  gap  width)  of  7855.  Furthermore,  transitional  values  of  7880  and 
7932  were  predicted  analytically  by  Vest  and  Arpaci  (1969)  and  Korpela 
et  aK  (1973),  respectively,  for  a  similar  vertical  slot  geometry. 

Lee  and  Korpela  (1983)  also  reported  a  wavenumber  of  2.80 
associated  with  the  six-cellular  flow  pattern  in  their  narrow  vertical 
slot.  They  defined  the  wavenumber  as  2tt/Z,  where  Z  represented  the 
center-point  distance  between  two  adjacent  cells.  When  this  study’s 
initial  seven-cellular  state  at  t  ^  35.0  was  compared  to  their  narrow 
vertical  slot  geometry  (aspect  ratio  =15),  the  wavenumber  was 
determined  to  be  approximately  2.60--a  good  correlation. 

However,  as  previously  mentioned,  the  cellular  flow  within  the 
vertical  slot  geometry  was  steady,  whereas  in  this  study,  the 
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multicellular  flow  at  G  =  5.2  seemed  to  vacillate  periodically  about 
the  seven-cell  state.  That  is,  upon  increasing  t,  a  cyclic  pattern 
resulted:  an  additional  eighth  cell  formed  on  top  of  the  80°  arc, 
reverted  back  to  seven  cells,  then  an  additional  eighth  cell  formed 
on  the  bottom  of  the  80°  arc,  and  again,  reverted  back  to  seven  cells, 
forming  an  8-7-8-7  cellular  pattern  (see  Figure  6.17).  In  the 
seven-cell  formation,  the  stronger  cells  appeared  near  =  90°  (see 

Figures  6.16  and  6.18).  The  remaining  cells  successively  decreased  in 
strength  (in  approximate  proportion  to  sin  4>)  as  they  proceeded  away 


stream  function  at  r  =  .5 /i>  =  90°  for  the  half 


from  the  90°  point. 


At  the  start  of  the  cellular  cycle,  the  seven  cells  weakened  (the 
maximum  stream  function  decreased)  as  the  eighth  cell  formed  on  top 
(see  Figure  6.19).  Then,  the  eight-cell  pattern  grew  in  strength  until 
the  two  larger  cells  near  90°  merged  into  one  with  the  maximum  stream 
function  shifting  upward  to  92°,  thus  forming  an  upward-shifted  seven 
cellular  flow  pattern.  Similarly,  these  seven  cells  then  weakened  as 
the  eighth  cell  formed  on  the  bottom  portion  of  the  chain  of  cells 
(see  Figure  6.20).  Again,  the  eight-cell  pattern  grew  in  strength  until 
the  two  larger  mils  merged  into  one  with  the  maximum  stream  function 
shifting  downwards  now  to  87'.  thus  forming  a  downward-shifted 
seven-cel lular  flow  pattern  completing  the  8-7-S-7  cycle.  A  cyclic 
henav'or  of  the  average  inner  shear-stress  also  accompanied  this 
■  tianging  el  lular  pattern.  Table  6.2  outlines  a  typical  cycle. 

Note  that  these  re  1  a 1 1 onsh i ps  would  repeat  f or  subsequent  cycles.; 


Table  b . 2 . 

characteristics  of  the  H- 

7-8-7-  cel  1 ular  cycle 

Ce  1  1 

t 

Average 

Maximum 

Maximum 

structure 

i  nner 

stream  function 

stream 

shear-stress 

location  function  value 

8- top 

35.0 

.23722 

87° 

.4506 

7-top 

280.0 

.23644 

9  2° 

.4815 

8-bottom 

350.0 

.237247 

82" 

.4501 

7-bottom 

680.0 

.23627 

O 

r- 

co 

.4824 

This  table,  along  with  the  particular  nature  of  the  cellular  activity, 
seems  to  imply  a  potential  energy  exchange  in  the  8-7-8-7  cellular  cycle. 
That  is,  the  most  potential  energy  seems  to  be  stored  within  the  seven¬ 
cell  formation,  corresponding  to  a  minimum  average  shear-stress  and  a 
larger  maximum  stream  function.  It  then  appears  that  some  of  the 
potential  energy  is  released  and  is  converted  to  shear  energy  which  aids 
in  the  creation  of  an  eighth-cell.  This  then  causes  a  rise  in  the 
average  inner  shear-stress  and  a  corresponding  decrease  in  the  maximum 
stream  function. 

The  strength  of  the  stream  function  at  ^  =  90°/r  =  .5  behaved  in 
a  square-wave-like  time-periodic  fashion,  as  clearly  illustrated  in 
Figure  6.17.  The  period  of  this  cyclic  pattern  was  approximately  800 
time  (t)  units. 
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The  fact  that  the  cellular  behavior  at  G  =  5.2  is  unsteady  can  be 
supported  by  the  work  of  Thomas  and  Vahl  Davis  (1970).  They  observed 
a  multicellular  flow,  similar  to  that  reported  by  Elder  (1965)  and 
Lee  and  Korpela  (1983),  in  narrow  vertical  cylindrical  annuli  (see 
Section  2.6).  Even  though  they  were  only  interested  in  predicting 
the  onset  of  this  multicellular  flow,  they  did  report  that  the  flow 
was  unsteady  and  attributed  this  unsteadiness  to  curvature  effects. 

Figure  6.21  displays  the  vorticity  distribution  at  the  inner 
cylinder  wall  for  both  unicellular  and  multicellular  flows.  The 
unicellular  steady-state  result,  for  G  =  5.1,  behaved  in  a  sinusoidal 
manner  with  the  maximum  vorticity  occurring  near  =  90°.  Then,  due  to 
the  multicellular  instability  at  G  =  5.2,  the  vorticity  distribution 


1-cell  (G=5.1):  steady 
7-cells  (G=5.2) :  t=1060.0 


Angular  variation  of  the  inner-wall  vorticity  for  various  cellular 
patterns 


was  clearly  disturbed  and  behaved  in  an  intense  oscillatory  manner  near 
the  vertical  portion  of  the  half-annulus.  This  extremely  wavy 
vorticity  pattern  shifted  only  slightly  as  the  number  of  .ells  ma>.-  v*- 
transition  from  ei ght-cel  1  s-top ,  to  seven-cells,  to  e'.;-'*:-  e‘ •• 


The  seven-cell  solution  experienced  the  strongest  variation  m  Wci 
vorticity.  For  individual  comparison,  the  four  curves  .  '  t  • «  .: 


Figure  6.21  are  depicted  separately  in  Figures  6.2.'  tnr  /  -  . 

To  examine  the  change  in  cellular  oenavior  w’t^  ;  ' 

5.5  was  calculated  using  the  6  =  5.2  seven-cellular  •••■*•.■ 
at  t  =  35.0.  Figure  6.26a  displays  the  time  vdrsi*:  -  •  •>. 

function  at  ^  =  90°/r  -  .5,  using  a  constant  *  •,  me  ,  t  to . 
to  t  =  400,  a  hint  of  double-periodicity  seems  t.  :<» 
afterwards,  the  flow  seems  to  become  aoer ' oj ’ . 
of  chaotic  motion.  The  multicellular  pattern  wa  -  >-  • 

G  =  5.2,  except  in  this  case,  the  cellular  ‘low  ■ me;  • 

between  nine,  ten  and  eleven  cells.  The  man 'mum  t  •  ,,  • 
occurring  near  'i>  =  90°,  was  greatest  for  the  nme-  o' 
least  for  the  eleven-cell  structure  A  typical  eleven-  . 
is  shown  in  Figure  6.26b.  Also,  in  general,  the  ell./  m 
G  =  5.5  was  similar  to  that  experienced  with  G  -  • 

The  preceding  results  show  that  a  mul  ticel  lular 
instability  is  indeed  probable  in  the  vertical  sections  .  * 
horizontal  annuli.  It  appears  that  this  instability  arises 
of  the  coupling  of  the  nonlinear  terms  with  the  buoyancy  term. 
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8-CELLS-BOTTOM 


Figure  6.25.  Angular  variation  of  the  8-cell -bottom  inner-wall  vorticity  for  G  -  5.2 
and  t  =  1210.0 
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-G  sin  *P,  in  the  vorticity  equation.  Considering  this  statement 

A 

together  with  the  fact  that  periodicity  is  lost  for  G  =  5.5,  one  might 
infer  that  the  effects  of  the  nonlinear  terms  become  more  and  more 

A 

pronounced  as  G  increases. 

Based  on  this  reasoning,  it  is  anticipated  that  for  some 

A 

intermediate  value  of  G  between  5.2  and  5.5,  the  possibility  of  an 
emerging  doubly-periodic  solution  is  likely,  thus  causing  the  system 
to  make  a  further  transition  on  its  ordered  path  to  chaos. 

6.2.2.  Pr  •»  0  perturbative  solution 

To  verify  the  Pr  ■>  0  perturbative  expansion  solution,  the 
analytical  results  will  be  compared  to  related  numerical  results  for 

A 

the  pretransitional  case  G  =  5.1.  As  discussed  in  Section  4.6,  only 
F1  and  (as  given  in  Eqs.  (4.93)  and  (4.94),  respectively)  appear  to 
be  significant  prior  to  the  onset  of  multicellular  flow.  This  claim 
will  now  be  tested  in  the  following  comparisons. 

An  analytical  expression  for  the  inner  wall  vorticity  was  found 
by  evaluating  at  r  =  0  for  G  =  5.1: 

W.  =  .425  sin  *  . 

inner 

A 

A  comparison  to  the  numerical  result  at  G  =  5.1  is  given  in  Figure 
6.27.  The  largest  deviations  occurred  near  <P  =  45  and  135°, 
corresponding  to  relative  errors  of  nearly  five  percent.  At  the  point 
of  maximum  vorticity,  'P  =  90°,  the  error  was  only  about  .1  percent. 
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Figure  6.27.  Angular  variation  of  the  steady-state  inner-wall  vorticity  for 


Similarly,  an  analytical  expression  for  the  vorticity  at  the 


horizontal  center-line  was  obtained  by  evaluating  W-|  at  'P  =  90°  and 
G  =  5.1: 


Wgoo  =  2.55  (r2  -  r  +  I)  . 

This  result  was  compared  to  the  numerical  result  in  Figure  6.28.  The 
vorticity  was  positive  near  the  walls  of  the  annulus,  .81  <  r  <  .19, 
signifying  the  boundary-layer  region.  It  became  negative  within  the 
center  portion  of  the  annulus,  representing  the  inviscid-core  region. 

A  maximum  error  of  approximately  3.5  percent  occurred  near  the 
switching  point,  where  the  vorticity  changed  from  positive  to  negative 
The  maximum  vorticity  occurred  at  the  walls,  r  =  0  and  1,  and  took  on 
a  value  of  approximately  .425.  For  r  =  .5,  the  vorticity  reached  a 
minimum  at  -.2125. 

The  analytical  expression  for  stream  function  at  =  90°  and 

A 

G  =  5.1  is  given  by: 


Fg0o  =  5.527125  r2  (r  -  l)2  . 

Figure  6.29  displays  its  comparison  to  the  numerical  result.  A 
maximum  error  of  about  one  percent  occurred  at  r  =  .5,  where  the 
stream  function  reached  its  maximum  value  of  approximately  .346. 

The  above  results  indicate  that  the  pretransi tional  flow  field 
(G  <  5.1)  can  be  adequately  represented  by  the  first-term  expansions 


on  of  the  stream  function  at  90°  for 


(W^  and  )  of  the  perturbation  solution,  which  do  not  account  for 
nonlinear  effects.  That  is,  the  problem  appears  to  be  linear  and  steady 

A 

up  to  G  =  5.1.  This  suggests  that  the  full  nonlinear  influence  only 

becomes  evident  just  prior  to  and  beyond  the  point  of  instability.  For 

this  case,  the  nonlinear  effects  appear  to  come  into  play  for  G  >  5.2. 

Figure  6.30  was  generated  in  order  to  show  how  the  Pr  +  0 

perturbative  solution  breaks  down  as  the  Prandtl  number  increases. 

1/4 

Using  the  relationship  G  =  G  Pr  ,  two  numerical  results  corresponding 
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to  Pr  =  .7  and  2.0  were  compared  to  the  analytical  solution  for  G  =  4.0. 
Surprisingly,  for  Prandtl  numbers  as  large  as  2.0,  the  comparison 
between  the  analytical  and  numerical  results  was  fairly  close.  For 
Pr  =  .7  and  2.0,  respective  maximum  errors  of  3.5  and  10.5  percent 
occurred  near  ^  =  140°.  Thus,  for  air,  the  use  of  Pr  +  0  expansions 
should  provide  for  relatively  good  estimates  of  the  local  flow 
variables,  at  least  for  the  pretransitional  cases. 

The  finite-Prandtl  number  perturbative  solution  of  Eq.  (4.29)  can 
also  be  used  to  obtain  analytical  results  for  small  Prandtl  number  cases. 
Using  the  three-term  perturbative  solution  (F^,  F?  and  F^),  all  of  the 
important  physics  can  be  taken  into  account,  including  the  first  effects 
of  nonlinearity  through  F^  and  F^.  For  liquid  mercury,  Pr  =  .02,  a 
six-cellular  counter-rotating  flow  was  encountered  at  G  =  3.5.  In  the 
complete  annulus,  the  large  kidney-shaped  cells  maintained  themselves 
near  =  90  and  270°,  while  two  smaller  counter-rotating  cells  formed 
first  at  the  bottom,  then  at  the  top  portions  of  the  annulus  -  V  =  0  and 
180°.  This  multicellular  flow  pattern  was  more  fully  developed  for 


gular  varia 


G  =  4.0,  and  is  depicted  in  Figure  6.31.  This  particular  flow  pattern 
was  very  similar  to  that  observed  by  Mack  and  Bishop  (1968)  and  Huetz 
and  Petit  (1974)  for  a  larger  gap  width,  G  =  1,  and  Pr  =  .02.  They 
graphically  depicted  the  multicellular  behavior  at  Ra  =  300,  but 
mentioned  that  it  initially  occurred  (analytically)  at  a  lower  value  of 
Rayleigh  number.  In  addition,  the  presence  and  shape  of  the  secondary 
cells  for  G  =  1  and  Pr  =  .02  were  confirmed  by  the  numerical  experiments 
of  Charrier-Mojtabi  et  al_.  (1979),  where  they  reported  a  similar  type  of 
multicellular  flow  that  began  to  develop  near  Ra  =  300. 

~  i/4  ~ 

Using  the  relation  G  =  Ra  G,  G  =  4.0  translates  into  a  Rayleigh 
number  of  256  for  G  =  1.0,  thus  reaffirming  the  results  described  above. 
Note  also  that  Huetz  and  Petit  (1974)  claimed  that  the  bottom  cell 
developed  first.  However,  when  the  finite-Prandtl  number  boundary-layer 
equations  were  solved  numerically  for  Pr  =  .02  (thus  incorporating  the 
full  effects  of  nonlinearity),  a  multicellular  flow  field  resembling 
that  described  with  the  Pr  -*■  0  numerical  solution  resulted  in  the 
vertical  section  of  the  half-annulus  at  G  =  1.92  (G  =  5.106). 

From  these  results,  one  might  conclude  that  the  effects  of  the 
boundary-layer  are  less  significant  in  the  larger  size  gaps,  thus 
creating  a  weakly  nonlinear  instability  as  originally  reported  by 
Mack  and  Bishop  (1968).  But  for  the  narrow-gap  widths,  the  boundary- 
layer  effects  are  much  more  dominant,  thus  creating  stronger  nonlinear 
interactions  which  trigger  the  hydrodynamic  type  of  instability  as 
described  in  this  study  for  the  small -gap/smal 1-Prandtl  number  limiting 
conditions.  These  results  are  discussed  further  in  Appendix  B. 


Figure  6.31.  Streamlines  of  the  steady  six-cellular  perturbation 
solution  corresponding  to  Pr  =  .02  and 
G  =  4.0  (diameter  ratio  =  2.0) 


6.2.3.  Small -Prandtl  number  stability  curves 

In  order  to  obtain  a  more  universal  relationship  for  predicting 
multicellular  flow  instability  with  regard  to  small -Prandtl  number 
fluids  in  the  narrow  horizontal  annuli,  the  following  derivation  was 
considered. 

It  was  shown  in  Chapter  4  that 


G  =  Ra1/4  G  =  Pr1/4  G 


(6.1) 


but  since  Ra  =  Gr  •  Pr,  Eq.  (6.1)  becomes  (after  some  manipulation): 
d 

Gra  =  (|)4  (6.2) 


where  Gr  is  the  Grashof  number  based  on  the  inner  cylinder  radius, 
d 

Then,  using  the  relation 


one  obtains: 


(6.3) 


where  Gr,  now  represents  the  Grashof  number  based  on  the  gap  width. 

D”d 

For  Pr  =  0,  the  transition  to  mul ticel lular  flow  occurred  at  G  =  5.2. 


Using  this  result  transforms  Eq.  (6.3)  into  a  universal  relation  which 


197 


is  inversely  proportional  to  the  gap  number,  G: 


(6.4a) 


This  relationship  is  plotted  for  various  values  of  G  in  Figure  6.32. 

Note  that  the  transitional  value  of  Gr^_a  =  8000,  as  predicted  in  this 
study  and  by  Lee  and  Korpela  (1983),  corresponds  to  a  gap  number  of 
approximately  G  =  .10.  Also,  the  shape  of  the  stability  curve  seems  to 
make  sense,  since  for  the  vertical  slot  studies  it  was  reported  by  Lee 
and  Korpela  (1983)  and  Elder  (1965)  that  as  the  gap  width  increased,  the 
transition  to  turbulence  occurred  sooner.  Similarly,  in  this  study,  it 
was  shown  for  very  small-gap  widths  (for  air)  that  as  the  gap  number 
decreased,  the  transitional  Rayleigh  number  for  multicellular  flow 
increased.  The  above  statements  indicate  that  as  gap  width  decreases, 

A 

the  onset  of  instability  is  suppressed.  Thus,  the  stress  parameter  G 
or  Rayleigh  number,  must  increase  accordingly  in  order  to  trigger  the 
respective  hydrodynamic  or  thermal  instability.  Qualitatively,  the  shape 
of  the  stability  curve  in  Figure  6.12  for  Pr  =  .706  (in  the  small -gap 
limit)  is  similar  to  that  in  Figure  6.32  for  Pr  =  0.  Also,  for  thermal 
instability,  the  following  relation  holds: 

r4 

GVa  =  G^TF 


where  for  air  in  the  small-gap  limit,  Pr  =  .706  and  Gyp  -  4.0,  hence: 


Gr 


b-a 
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G 


air 


(6.4b) 


Prandtl  number  stability  curve:  Grashof  number 
,  Gr.  ,)  as  a  function  of  gap  number 


This  relationship  is  very  similar  to  that  for  hydrodynamic  instability 
(Pr  =  0),  except  the  proportionality  constant  differs  by  a  factor  of 
two. 

Taking  the  above  approach  a  step  further,  one  can  also  determine  a 

stability  curve  based  on  fini te-Prandtl  number  by  using  the  relation 
~  1/4 

G  =  5.2  Pr  ,  where  G  =  5.2  again  represents  the  transitional  value  for 
Pr  =  0.  The  new  stability  curve  is  shown  in  Figure  6.33.  These 
analytical  predictions  were  compared  against  numerically  predicted 
transitional  values  of  G  (obtained  from  the  finite-Prandtl  number 
boundary- layer  equations)  for  four  different  smal 1 -Prandtl  number  fluids 
(see  Table  6.3). 

Table  6.3  indicates  that  for  small-Prandtl  numbers,  the 
~  1  /4 

relationship  G  =  5.2  Pr  '  appears  to  be  valid  for  predicting  the 
transitional  G  associated  with  hydrodynamic  instability  in  the  vertical 
sections  of  narrow  annuli.  For  Pr  =  .02,  an  initial  seven-cell  unsteady 
flow  field,  resembling  that  of  Pr  =  0,  was  obtained  at  G  =  1.92. 

Table  6.3.  Numerical  predictions  of  G  transitional  as  a  function  of 
Prandtl  number 


Prandtl  number 


G  (transitional) 


G(=  G-Pr"1/4) 


e-Prandtl  number  stability  curve 
function  of  Prandtl  number 


However,  for  Pr  =  .35  (at  G  =  3.99),  an  interesting  change  occurre d-- 
only  a  six-cell  pattern  resulted  (see  Figure  6.34)  and  the  solution 
achieved  a  steady-state  condition  similar  to  that  reported  for  the 
cellular  instability  between  narrow  vertical  slots.  Also,  at  Pr  =  .35, 
the  isotherms  were  no  longer  perfectly  concentric  as  with  the  Pr  =  0 
solution.  Instead,  they  responded  quite  closely  with  the  fluid  motion 
(see  Figure  6.35). 

From  these  results,  it  appears  that  the  same  type  of  hydrodynamic 
instability  as  for  Pr  =  0  extends  to  at  least  Pr  =  .35.  But,  the  nature 
of  the  multicellular  flow  has  changed  from  seven-cells  unsteady  for 
Pr  =  .02  to  six-cells  steady  with  Pr  =  .35.  Perhaps,  though,  if  one 
stresses  the  Pr  =  .35  case  past  G  =  3.99,  an  unsteady  cellular  motion 
may  ensue.  Note  that  for  Pr  =  .2  and  .3,  the  cellular  flow  field  was 
not  resolved;  only  the  point  of  instability  was  determined. 

The  two  stability  curves  described  above  should  provide  for  simple 
and  convenient  comparisons  to  future  related  studies  of  narrow 
horizontal  annuli  (especially  experimental  ones).  Note  that  the  curve 
in  Figure  6.33  could  easily  be  extended  for  comparison  to  data  in  the 
larger  Prandtl  number  range,  thus,  eventually  establishing  an 
acceptance  band  on  this  asymptotically  generated  stability  curve. 


Figure  6.34.  Streaml 
for  Pr 
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7.  CONCLUSIONS 

Numerical  solutions  of  the  2-D  Navier-Stokes  equations  yielded 

steady  laminar  multicellular  flow  for  air  near  the  top  portions  of 

narrow  horizontal  isothermal  concentric  cylinders.  It  was  shown  that 

hysteresis  behavior  exists  for  G  =  .100  and  .200,  with  regard  to  sudden 

changes  in  the  mean  Nusselt  number  and  average  shear-stress  as  the 

Rayleigh  number  was  slowly  increased  and  then  decreased  past  the 

multicellular  transition  point.  Also,  for  G  less  than  approximately 

.175,  the  results  for  air  indicate  that  the  transitional  Rayleigh  number 

increases  as  the  gap  number  decreases.  The  minimum  transitional  Rayleigh 

number  occurred  at  Ra.  ,  =  2400  for  G  =  .175. 

b-a 

In  performing  this  multicellular  flow  study,  first-order  upwind 
differencing  of  the  nonlinear  convective  terms  proved  inadequate.  In 
relation  to  the  pretransitional  flow  field  cases,  the  first-order 
method  was  in  good  agreement  with  the  results  obtained  from  the 
second-order  central -differenced  numerical  scheme.  However,  for  the 
higher  Rayleigh  number  flow  regimes,  the  artificial  viscosity 
associated  with  the  first-order  differencing  of  the  nonlinear  convective 
terms  suppressed  the  development  of  secondary  counter-rotating  cells. 
Thus,  it  appears  that  at  least  second-order  accuracy  is  needed  to 
properly  capture  multicellular  flow  behavior.  Moreover,  multicellular 
flow  resolution  is  also  strongly  dependent  upon  mesh  size.  That  is, 
coarse  meshes  frequently  failed  to  resolve  the  multicells,  while 
finer  meshes  would  sometimes  result  in  asymmetric  cellular  patterns 


when  not  enough  nodes  were  concentrated  in  the  thermal  plume  region 
for  air. 

A  high  Rayleigh  number/small-gap  asymptotic  expansion  theory  was 
constructed  which  simplified  the  2-D  Navier-Stokes  equations  into 
Cartesian-like  boundary-layer  equations.  Analytical  perturbative 
solutions  to  these  equations  were  obtained  and  the  results  compared 
favorably  with  the  pretransitional  numerical  data  generated  in  the  2-D 
Navier-Stokes  hysteresis  study  of  this  thesis.  For  the  limit  as 
G  -*•  0,  the  boundary-layer  equations  reduced  to  viscous-dominated 
Stokes-flow  equations,  as  shown  in  Chapter  4.  This  implies  that  for 
either  Ra  +  0  or  G  +  0 ,  a  conductive-dominated  flow  results. 

In  addition,  simplified  limiting  equations  for  Pr  •>  0  and  Pr  ■+  °° 
were  obtained  from  the  finite-Prandtl  number  boundary-layer  equations. 
In  the  Pr  -*■  °°  limit,  only  the  nonlinear  terms  in  the  energy  equation 
remained  and  the  vorticity  equation  reduced  to  a  Stokes-flow  equation, 
suggesting  the  possibility  of  a  Benard-type  thermal  instability,  as 
seen  with  air  near  the  top  portion  of  narrow  annuli.  In  direct 
contrast,  the  energy  equation  decoupled  from  the  vorticity  equation 
and  reduced  to  simply  T  =  1-r  for  the  limit  as  Pr  0.  Thus,  only  the 
nonlinear  terms  in  the  vorticity  equation  came  into  play,  suggesting 
a  source  of  hydrodynamic  instability. 

When  the  Pr  -+  0  equations  were  solved  numerically,  an  unsteady 
time-periodic  multicellular  instability  developed  (at  G  =  5.2)  in  the 
vertical  portions  of  a  narrow  horizontal  annulus.  Both  the  wavenumber 
and  the  transitional  Grashof  number  (based  on  gap  width)  of  the 


206a 


initial  seven-cellular  instability  compared  favorably  with  the  steady 
six-cellular  instability  experienced  in  a  narrow  vertical  slot  of 
aspect  ratio  equal  to  15,  as  studied  by  Vest  and  Arpaci  (1969)  and 
Lee  and  Korpela  (1983). 

Upon  increasing  time,  cyclic  changes  in  the  multicellular  flow 

/V 

pattern  resulted  for  G  =  5.2,  forming  an  8-7-8-7  repetitive  cellular 

A 

behavior.  However,  this  periodic  flow  pattern  was  lost  when  G  was 
increased  from  5.2  to  5.5,  resulting  in  a  seemingly  chaotic  cellular 
behavior  between  nine,  ten  and  eleven  cells. 

Also,  the  Pr  -*■  0  analytical  pretransitional  solutions  were  in 

/V 

good  agreement  with  the  numerical  results  up  to  G  =  5.1.  Based  on  this 
comparison,  the  Pr  -*-  0  problem  appeared  linear  and  steady  prior  to  the 
point  of  instability. 

Moreover,  a  universal  relationship  for  small-Prandtl  number 
fluids, 


was  suggested  as  a  simple  guide  for  predicting  the  onset  of 
hydrodynamic  instability  between  narrow  vertical  slot-like  geometries. 

Figure  6.13  shows  that  for  air  (Pr  =  .706)  in  very  narrow-gap 
widths  (G  <  .075),  thermal  instability  sets  in  near  G  =  3.90.  But  for 

/V 

Pr  =  0,  hydrodynamic  instability  sets  in  at  G  =  5.2,  which  corresponds 
to  a  G  of  approximately  4.8  for  air  when  using  the  analytical  relation 
G  =  5.2  Pr^.  This  difference  in  G  indicates  that  the  mechanism  which 


triggers  the  instability  for  air  is  different  than  that  for  Pr  =  0  (as 

previously  discussed),  and  is  probably  largely  due  to  the  nonlinear 

terms  in  the  energy  equation  being  coupled  with  the  buoyancy  term 
3T 

cos  P  in  the  2-D  Navier-Stokes  equations  (see  Appendix  B). 
Considering  the  above,  competing  effects  of  both  thermal  and 
hydrodynamic  instability  may  result  for  Prandtl  numbers  between  zero 
and  .706. 

In  order  to  more  fully  support  and  study  the  nature  of  the 
unsteady  multicellular  instability  described  in  this  research  effort, 
an  experimental  analysis  of  a  narrow  horizontal  annulus  filled  with 
liquid  mercury  (Pr  =  .02)  is  recommended  as  a  possible  follow-up 

A 

project.  Also,  numerical  solutions  relating  to  G  =  5.3  or  5.4  should 
be  carried  out  to  determine  if  a  doubly-periodic  cellular  behavior  does 

A 

indeed  precede  the  chaotic  state  noted  at  G  =  5.5,  as  suggested  in 
the  'ordered  path  to  chaos'  scenario  described  by  Ruelle  and  Takens 
0971  ). 

Ay 

In  addition,  the  energy  spectrum  of  the  fluid  motion  at  G  =  5.2 
and  5.5  could  be  calculated  in  order  to  reveal  any  type  of  spectral 
broadening  (or  fundamental  frequencies)  associated  with  flow  changes 
from  time-periodic  to  chaotic-like  behavior. 

Lastly,  the  Pr  -+■  00  boundary-layer  equations  should  be  solved 
numerically  to  determine  the  nature  of  any  cellular  instability  that 
might  develop.  In  the  same  spirit,  the  limiting  2-D  Navier-Stokes 
equations  derived  in  Appendix  B  should  also  be  further  investigated 
numerically  for  different  size  gap  widths. 
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In  closing,  since  natural  convective  flow  instabilities  generally 
develop  slowly,  especially  in  narrow  horizontal  annuli  where  the 
geometry  itself  limits  the  number  of  degrees  of  freedom,  a  more 
thorough  investigation  of  laminar  flow  transition  can  be  undertaken. 
Also,  it  has  been  conjectured  that  the  transition  to  chaos  in  a  system 
with  many  degrees  of  freedom  (such  as  a  fluid  flow  at  large  Reynolds 
numbers)  is  not  qualitatively  different  from  the  transition  to  chaos 
in  a  system  with  a  few  degrees  of  freedom  (Jensen,  1987).  Thus,  it  is 
hoped  that  further  study  of  simple  geometries,  such  as  described  in 
this  thesis,  will  lead  to  a  deeper  understanding  of  the  nonlinear 
processes  involved  in  the  mysterious,  yet  ever-present  dynamics  of 
fluid  flow  instability  phenomena. 
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10.  APPENDIX  A:  ANALYTICAL  COEFFICIENTS 


The  perturbative  coefficients  corresponding  to  the  analytical 
solution  of  the  finite-Prandtl  number  boundary-layer  equations, 
(4.39)  -  (4.46),  are  listed  below.  Note  again  that  all  of  the 
coefficients  (A2  through  H4)  are  functions  of  r  only. 
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11.  APPENDIX  B:  LIMITING  CONDITIONS  OF  THE  2-D 
NAVIER-STOKES  EQUATIONS 

Since  the  simplified  boundary-layer  equations  of  Chapter  4  (in 
the  limits  of  Pr  -*■  0  and  Pr  ®)  applied  only  to  the  narrow-gap/high 
Rayleigh  number  regime,  the  effects  of  small  and  large  Prandtl  number 
on  arbitrary  gap  size  and  Rayleigh  number  could  not  be  explored.  To 
consider  these  effects,  the  2-D  Navier-Stokes  equations  of  Chapter  3 
must  be  employed.  This  appendage  will  derive  these  limiting  equations, 
and  will  explain  the  important  features  and  differences  of  both 
approaches. 


T(0,  'P)  =  1  (11.2a) 

T ( 1 ,  'f')  =  0  (11.2b) 


U  (r,  *  =  0,  IT)  =  0 


(11.2c) 
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Given  Eqs.  (11.2a-c),  the  solution  to  Eq.  (11.1)  is 


T  = 


£n  ( 1  +  G  •  r ) 
£n  (1  +  G) 


(11.3) 


which  is  simply  the  steady-state  conduction  solution  for  the  concentric 
cylindrical  annulus.  Because  Eq.  (11.3)  is  independent  of  one  gets: 


(11.4) 

and 

9T  -  "2  n 

3r  '  (1  +G1  r  )£n(l  +  G)  * 

Hence,  for  the  limit  as  Pr  0,  the  buoyancy  term,  -cos  <p  in 
Eq.  (3.1b),  reduces  to  zero,  whereas  sin  'P  becomes: 


sin  *  |I  - 


IT 


-G  sin  4> 
G  ■  rj  £n  ( 1 


■GT 


(11.6) 


Thus,  the  energy  equation  has  again  decoupled  itself  from  the  vorticity 
equation. 

Now,  in  the  limit  as  G  -*•  °°,  Eq.  (11.6)  vanishes  (or  approaches 
zero  from  the  negative  side).  But  as  G  0,  Eq.  (11.6)  simplifies  to: 

=  -sin  'P  .  (11.7) 

G->0 


sin  * 

3r 


,*  V  .  -  •  *  -  v  ‘O  *  V 
■a  ^  ^ 
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9T 

Therefore,  sin  *P  ^  achieves  its  maximum  absolute  value  at  ^  =  90  and 


270°  for  G  0,  thus  duplicating  the  Pr  -*■  0  condition  described  in 


Chapter  4  in  which  a  hydrodynamic  (or  shear-flow)  instability  is 


probable  for  the  narrow-size  gap  widths. 


From  the  above  results,  one  may  deduce  that  as  G  -*•  0,  the  shearing 


influence  of  the  inner  and  outer  boundary-layers  is  strong  enough  to 


cause  hydrodynamic  instability  to  originate  near  the  points  of  maximum 


velocity,  or  at  ^  =  90  and  270°.  In  contrast,  as  gap  size  increases. 


the  shearing  influence  of  the  boundary-layers  is  less  destabilizing 


since  the  buoyancy- induced  velocities  are  reduced.  However,  at  high 


enough  Grashof  numbers  within  the  inviscid  core,  a  thermal-like 


multicellular  instability  develops  near  the  upper  and  lower  portions  of 


the  annulus  (^  =  0  and  180°),  where  smaller  counter-rotating  cells  have 


been  reported  (Mack  and  Bishop,  1968;  Huetz  and  Petit,  1974;  Charrier- 


Mojtabi  et  al. ,  1979)  (also  see  Subsection  6.2.2  and  Figure  6.31). 


After  factoring  out  the  Prandtl  number  dependency  in  the  vorticity 


equation  of  (3.1b)  and  using  the  results  of  Eqs.  (11.4)  and  (11.5), 


the  governing  equations  for  the  Pr  ^  0  limit  become: 


Vorticity: 


p2  3w  ,  /  ,  1 \  / 3f  3w  9f  3W\ 
G  at  (r  G}  W  aF' 


2-?+  (r  + 1'"'  I? +  (r  +  s> 2 

3r2  S  3r  G  34,2 


,  r/n  \  r  -G  sin  3 

+  G(Gr)  {(1  +  G  •  r)  &n  ( 1  +  G)} 


(11.8a) 


Stream  Function: 


4+(r  +  -)-1|l+(r  +  l,-23!f,G2, 

3r  b  9r  G  3^ 


(11.8b 


with  boundary  conditions: 


f(0,  *)  =  f (i ,  n»)  =  o 


(11.8c 


2 

w(r  =  0,  1;  40  =  V 

6  3r  r=0,l 


(11. 8d 


Consistent  with  the  Pr  -*■  0  result  discussed  in  Chapter  4,  the  nonlinear 
terms  in  the  vorticity  (momentum)  equation  remain  and  the  energy 
equation  decouples  and  reduces  to  the  simple  conduction  solution. 
Importantly,  though,  the  equations  given  in  (11.8)  are  valid  for 
arbitrary  gap  and  Grashof  number,  thus  enabling  one  to  further 
investigate  the  various  effects  of  each  parameter. 


11.2.  The  Pr  ■*  °°  Equations 

For  the  boundary-layer  equations  of  Chapter  4,  in  the  limit  as 
Pr  -+  ",  specific  expansions  were  derived  in  order  to  retain  all  of 
the  important  physics.  Using  the  same  rationale  (see  Eq.  (4.102)), 
the  following  expansions  will  be  adopted  for  this  more  general  approach 
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«*.  «■  / 


w  =  Pr  W  +  0(Pr  ) 


f  =  Pr'1  F  +  0(Pr'2)  (11.9) 

T  =  T  +  0(Pr~] )  . 


Then,  in  the  Pr  -*■  00  limit,  Eqs.  (3.1  a-e )  simplify  to: 


Thermal  Energy: 
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Vorticity : 
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(11.10b) 


Stream  Function: 


-1  ~  -2 

3^F  .  ,  A  1\  9F  +  #r  .  lx  3 F 
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(11.10c) 


with  boundary  conditions: 


T(0,  40  =  1;  T(l, 

)  =  o 

(ll.lOd) 

F ( 0 ,  40  =  F(l,  40  = 

0 

(ll.lOe) 

W(r  =  0,  1;  40  = 

rF 

2 

(11.1  Of) 

r=0,l 

In  the  above  equations,  only  the  nonlinear  terms  in  the  thermal -energy 
equation  remain  and  the  Prandtl  number  dependency  has  again  vanished. 
However,  both  buoyancy  terms  are  retained  in  the  Stokes-like  vorticity 
equation,  and  Eqs.  (11.1 Oa-f )  are  also  valid  for  arbitrary  gap  and 
Rayleigh  number. 

As  Pr  -*  oo,  the  thermal -energy  diffuses  much  more  slowly  than  for 

9T 

Pr  -*•  0;  hence,  the  radial  temperature  gradient,  — ,  tends  toward  its 

o  r 

minimum  in  the  former  case  and  approaches  its  maximum  in  the  latter. 

This  implies  that  as  Pr  -*■  °°,  the  tendency  of  a  hydrodynamic  instability 

to  originate  near  the  vertical  portion  of  the  annulus  is  reduced, 

9T 

since  the  term  sin  is  rot  as  strong  as  in  the  Pr  -*■  0  case.  This 

dr 

statement  is  further  supported  by  the  results  of  Lee  and  Korpela  (1983). 
They  found  that  as  Prandtl  number  increased,  the  width  of  the 
vertical  slot  had  to  be  decreased  (or  aspect  ratio  increased)  in  order 
to  trigger  the  multicellular  instability  more  readily  seen  in  the 
smaller  Prandtl  number  fluids. 


Furthermore,  regardless  of  the  flow,  only  concentric  isotherms 


result  for  Pr  0.  However,  the  isotherms  are  significantly 

influenced  by  the  flow  field  for  Pr  -►  thus  causing  the  angular 

9T 

temperature  gradient  ^  to  be  of  much  more  importance.  This  behavior 

is  especially  true  for  the  thermal  plume  region  near  the  top  of  the 

inner  cylinder  =  180°),  whereas  near  the  bottom  of  the  annulus 

(4>  =  0°),  tt  is  minimal  for  flow  fields  in  this  larger  Prandtl  number 

range.  Then,  since  cos  4Ms  -1  at  4*  =  180°  and  zero  at  4>  =  90°, 

and  because  is  most  significant  in  the  vicinity  of  'P  =  180°,  the 

buoyancy  term  cos  tp  reaches  its  absolute  maximum  near  V  =  180°, 

where  sin  —  tends  to  zero, 
dr 

From  the  above  discussion,  the  following  conclusions  can  be 


drawn.  Due  to  the  nonlinear  terms  in  the  energy  equation  being  coupled 
with  the  buoyancy  term  cos  4>  a  thermal -type  multicellular  flow 
instability  near  the  top  of  narrow  annuli  may  occur  for  large  Prandtl 
number  fluids.  Whereas,  for  Pr  -*■  0,  a  hydrodynamic  instability  is 
probable  near  the  vertical  portions  of  narrow  annuli,  due  to  the 


nonlinear  terms  in  the  vorticity  equation  being  coupled  with  the 

piT 

buoyancy  term  sin  ^  — ,  as  previously  shown  in  Chapter  4. 


12.  APPENDIX  C:  E  VERSUS  E  '  FORMULATION 


As  discussed  in  Chapter  5,  a  correction  term  was  added  to  the 
first-order  upwind-differenced  representations  of  both  the  streamwise 
(angular)  and  transverse  (radial)  convective  terms.  This  was  done  in 
order  to  gain  second-order  central -differencing  accuracy  while  avoiding 
the  stability  problems  of  spurious  wiggles  at  high  Rayleigh  numbers, 
associated  with  ordinary  central -di fferencing  of  the  nonlinear  terms. 

For  the  2-D  Navier-Stokes  equations,  the  correction  term  (E)  in 
Eqs.  (5.6a-f)  was  evaluated  at  the  old-time  level  En,  as  opposed  to  the 
new-time  level  En  ,  for  a  truly  implicit  method.  This  particular 
change  seemed  to  further  aid  in  stabilizing  the  numerical  method  for 
highly  convective  flows,  and  also  significantly  reduced  the  total 
number  of  iterations  required  for  convergence  to  a  steady-state 
solution,  as  will  be  shown  in  Section  12.2. 

To  show  that  this  substitution  (En  for  En+^ )  does  not  affect  the 
formal  accuracy  of  the  numerical  method,  a  formal  truncation-error 
analysis  was  performed  on  the  governing  set  of  equations  as  written 
in  (5.4)  and  (5.6).  The  basic  procedure  and  final  results  are  given 
in  the  following  section. 

12.1.  En  Truncation  -  Error  Analysis 
Using  the  En  formulation,  the  general  finite-difference  equation 
for  4)  =  T,  w  or  f  can  be  written  as: 
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(12. 2e 


The  remaining  terms  in  En,  <}>",  and  are  represented  by  double 

Taylor-series  expansions;  for  example: 
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(12. 2f 
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and  similarly  for  ^  ^  and  4>^. 

Considering  the  stream  function  equation,  the  appropriate 
expansions  of  Eq.  (12.2)  are  substituted  into  Eq.  (12.1),  resulting 
in  the  following  simplified  expression: 


n+1  2 

+  <r  +  h  % 

n+1 

<f> 

rr 

+  2Xd>r 

0 

0 

n+1 

0 


+  (TE)f  =  0  (12.3) 

f 

J 

where  (TE)^  signifies  the  stream  function  truncation-error  and  is 
given  by: 

(TE)f  =  0[At(hf  -  hb),  hf  -  hb,  kf  -  kfa,  hf2,  hb2, 

At3,  kf3,  kb3]  . 

Note  that  for  the  stream  function  equation,  both  the  unsteady  and 
streamwise  convective  terms  are  absent,  as  shown  in  Eqs.  (5.4)  and  (5.5). 
Thus,  Eq.  (12.3)  recovers  the  proper  differential  equation,  and  (TE)f 
identifies  that  the  formal  accuracy  of  the  differencing  method  is  not 
altered  when  replacing  En+1  with  En  in  the  stream  function  equation. 

In  a  similar  manner,  the  energy  equation  simplifies  to: 

i 


m  *  S  * 


;* . 
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(TE)t  =  0[At(hf  -  hfa),  At(kf  -  kb),  At,  hf  -  hb,  kf  -  kb. 


Hence,  the  original  energy  equation  is  recovered  and  the  formal 
accuracy  of  the  method  is  not  sacrified.  Likewise,  the  vorticity 
truncation-error  is  of  the  same  order  as  the  energy  equation. 

To  compare  the  differences  of  the  two  formulations,  several  test 
cases  were  performed  and  are  discussed  in  the  next  section. 

12.2.  Numerical  Comparison  of  En  Versus  En+1 

For  the  test  cases  that  follow,  steady-state  results  were  obtained 

-4 

by  using  a  time-step  of  5  x  10  carried  out  to  t  -  .150.  Also,  the 
0  -  2tt  program  was  used  with  a  mesh  size  of  31  x  102  nodes.  For  the 
pretransitional  cases  compared  below,  the  gap  size  was  G  =  .200  and 


Table  12.1.  Pretransition  comparison  of  En  versus  En+^ 


Rayleigh 

number 

E" 

iterations 

^n+1 

iterations 

Convergence 

criteria 

100,000 

4585 

4585 

1  x  10'6 

250,000 

6126 

6550 

1  x  10'6 

350,000 

9520 

10590 

X 

o 

1 

cn 

the  initial  conditions  were  zero  (see  Table  12.1).  At  Ra  =  100,000, 
the  flow  field  was  conduction-dominated  and  the  number  of  iterations 
for  convergence  was  the  same  for  both  formulations.  But  when  the 
flow  became  convecti ve-dominated  at  Ra  =  350,000,  the  En  formulation 
improved  the  convergence  rate  by  approximately  10  percent.  For  these 
three  pretransitional  comparisons,  both  methods  produced  exactly  the 
same  average  Nusselt  numbers  and  shear-stresses. 

However,  checking  one  of  the  hysteresis  loop  calculations  for 
G  =  .200,  a  marked  improvement  could  be  seen  in  the  number  of 
iterations  saved  for  the  four-cellular  initial  condition  (Ra  =  900,000), 
used  to  start  the  Ra  =  350,000  test  case.  The  results  are  summarized 
in  Table  12.2.  In  this  case,  the  convergence  rate  of  the  En  method  was 
approximately  25  percent  faster  than  the  En+^  approach.  For  the  En+1 
method,  the  two  cells  began  to  bottle-neck  (local  minimum  stream 
functions  resulted)  near  ^  =  170  and  190°,  but  never  split  apart  to 
form  the  four-cellular  counter-rotating  pattern.  But,  when  the  En+1 
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Table  12.2.  Multicellular  flow  comparison  of  En  versus  En+"' 


Rayleigh 

Iteration 

Convergence 

Cellular 

number 

number 

criteria 

structure 

En 

350,000 

5850 

1  x  10‘6 

4-cells 

^n+1 

350,000 

7790 

1  x  10"6 

2-cells 

method  was  repeated  with  a  convergence  constraint  of  7.5  x  10”^,  the 

n  +  «| 

4-cell  pattern  was  achieved.  Therefore,  it  appears  that  as  the  En 
convergence  constraint  is  reduced,  its  calculations  compare  more 
favorably  to  those  of  En  for  multicellular  flows. 

Thus,  the  En  method  has  both  speed  and  accuracy  advantages  when 
multicellular  flow  is  experienced.  This  is  most  probably  due  to 
added  stability,  since  En  is  always  constant  when  iterating  on  flow 
variable  calculations  for  a  new  time-step.  Also,  a  relaxation  strategy 
is  not  necessary  when  using  the  En  formulation.  For  these  reasons,  the 
En  approach  was  used  in  the  2-D  Navier-Stokes  numerical  scheme. 


13.  APPENDIX  D:  HEA  TRANSFER  AND  SHEAR 
STRESS  RELATIONS 

13.1.  Dimensionless  Nusselt  Numbers 
The  local  heat  transfer  rates  along  the  heated  inner  cylinder  wall 
and  the  cooled  outer  cylinder  wall  are  calculated  from  first  principles. 
For  the  inner  and  outer  cylinders,  these  local  rates  are  given  by: 


where  a,  b  are  the  radii  of  the  inner  and  outer  cylinders, 
respectively.  Using  the  definition  of  Nusselt  number  based  on  the 
inner  cylinder  radius,  Eqs.  (13.1)  and  (13.2)  can  be  nondimensional ized 
to  yield  the  following  expressions: 


(13.3) 

(13.4) 


The  average  heat  transfer  rates  can  be  found  by  simply  integrating 


the  local  radial  temperature  gradients  along  the  inner  and  outer 
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(13.5) 


-1 

2ttG 


2tt 

/ 

0 


d  ^  . 


(13.6) 


Note  that  when  symmetry  is  assumed,  the  integrations  are  carried  out 
for  only  the  half-annulus,  or  from  0  to  tt.  Numerically,  the  integrations 
of  (13.5)  and  (13.6)  are  approximated  by  using  a  second-order  trapezoidal 
rule  with  variable  increments,  thus  preserving  the  formal  second-order 
accuracy  of  the  numerical  method.  Also  note  that  since  there  are  no 
internal  heat-generating  sources  (or  sinks),  the  inner  and  outer 
average  Nusselt  numbers  should  equal  each  other  upon  reaching  steady 
state,  i .e. ,  Nu^  =  NuQ. 


13.2.  Dimensionless  Shear  Stresses 
Also  from  first  principles,  the  local  shear  stresses  at  the  inner 
and  outer  cylinder  walls  are  defined  as: 
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(13.7) 

(13.8) 


where  y  is  the  dynamic  viscosity  coefficient  and  —  is  the  radial 

9r 

gradient  of  the  tangential  velocity  component.  Expressing  v  in  terms 


of  the  stream  function  and  then,  nondimensionalizing  yields  the 
following  relations  for  the  local  shear  stresses: 
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(13.9) 


(13.10) 


The  average  shear  stresses  can  also  be  obtained  by  integrating  the 
local  values  along  the  inner  and  outer  cylinders: 
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(13.11) 

(13.12) 


Again,  the  integrations  were  carried  out  numerically  by  using  a 
second-order  trapezoidal  rule  with  variable  increments. 


13.3.  Heat  Transfer  and  Shear  Stress  Finite-Difference 

Expressions 

In  order  to  numerically  approximate  the  local  and  average  Nusselt 
numbers  as  given  in  Eqs.  (13.3)  -  (13.6),  finite-difference  expressions 
for  the  radial  temperature  gradient  on  the  inner  and  outer  cylinder 
walls  must  first  be  obtained.  To  do  this,  a  second-order  polynomial 
curve  fit  was  employed: 


Using  the  inner  cylinder  boundary  condition,  T?+]  =  1,  the  following 

'  >J 

equations  can  be  written: 


a^  =  1  for  all  j 
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where  Ar-|  =  0  and  Ari  s  r.  -  r^_1  for  i  =  2,  3.  Note  that  i  =  1 
corresponds  to  the  inner  cylinder  wall  at  r  =  0  and  i  =  2,  3  represent 
consecutive  nodal  points  in  the  positive  radial  direction.  Then,  since 


II 

dr 
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l,j 


one  obtains,  after  some  algebra,  the  following  expression  for  the 
inner  local  Nusselt  number: 
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where 
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T3+]  -  1  +  Ar  -  T2+]) 

n+1  _  _ _ 2 _ 

Cj  "  Ar 

<4r2  *  Ar3>  <’  -  A?2hr3> 


Similarly,  using  the  outer  wall  boundary  condition,  t!"!^  .  =  0,  the 

NK ,  J 

expression  for  the  outer  local  Nusselt  number  becomes: 


Nu0  = 


Tn+1  ArNR  +  ArNR1  n+1 
NR1 ,j  1  ArNR  >  "  NR2,j 

wArNR  +  ArNRl  ,, 
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I 

where 


ArNR  ‘  rNR  '  rNRl 
ArNRl  =  rNRl  '  rNR2 


(13.16) 


and  i  =  NR  corresponds  to  the  outer  cylinder  wall  at  r  =  1,  while 

i  =  NR1 ,  NR2  represent  consecutive  nodal  points  in  the  negative  radial 

direction. 

In  the  same  fashion,  a  similar  polynomial  fit  as  in  Eq.  (13.13) 
can  also  be  used  for  the  stream  function  in  order  to  obtain 
expressions  for  the  inner  and  outer  local  shear  stresses.  Employing 
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the  boundary  conditions 


fn+l  =  fn+l  =  0 
TU  NRJ  U  ’ 


Eqs.  (13.9)  and  (13.10)  can  be  cast  into  the  following  algebraic 
representations: 
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(13.18) 


The  above  relations  are  valid  for  the  2-D  Navier-Stokes  equations 
as  used  in  this  study.  However,  these  relations  can  also  be  used  with 
the  finite  and  zero-Prandtl  number  boundary-layer  equations  by  simply 
substituting  the  appropriate  related  variables  for  temperature,  stream 
function  and  gap  number.  Also,  since  a  second-degree  polynomial  curve 
fit  was  used,  the  Nusselt  number  expressions  are  typically  of  second- 
order  accuracy  for  uniform  increments,  while  the  shear  stress 
expressions  are  of  first-order  accuracy. 
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14.  APPENDIX  E:  TRUNCATION-ERROR  STUDY 


This  appendix  investigates  the  effect  of  truncation  error  on  the 
numerical  solution  to  the  2-D  Navier-Stokes  equations.  Work  begins 
with  the  comparison  of  first-order  upwind  differencing  versus  second- 
order  central  differencing  of  the  nonlinear  convective  terms  in  the 
2-D  Navier-Stokes  equations.  Then,  a  mesh  resolution  study  for  three 
different  size  gap  widths  is  conducted  to  gain  an  idea  of  the  uncertainty 
involved  in  the  numerically  predicted  transitional  Rayleigh  numbers  of 
Section  6.1.  Lastly,  two  correlations  for  the  mean  Nusselt  number  are 
compared  to  related  numerical  results  with  respect  to  an  annulus  of 
diameter  ratio  1.200. 

14.1.  Upwind  Versus  Central  -  Differencing 
In  Figure  14.1,  the  two  different  schemes  are  compared  for  a 
pretransitional  Rayleigh  number  of  300,000.  The  complete  annulus  was 
considered  with  G  =  .200  and  a  mesh  size  of  31  x  102  nodes.  The 
upwind  scheme  was  employed  by  setting  the  correction  term  (En)  to  zero 
in  Eq.  (5.6g).  The  first-order  upwind  method  compared  rather  well  to 
the  central -differencing  scheme  in  this  case,  with  maximum  relative 
errors  of  2-3  percent  occurring  near  ^  =  180°  and  V  =  0°,  respectively. 
Hence,  for  the  pretransitional  cases,  the  first-order  upwind- 
differencing  method  should  yield  fairly  accurate  numerical  results. 

However,  for  Ra  =  1,000,000  {well  beyond  the  point  of  multicellular 
transition  as  predicted  with  the  central-differencing  scheme),  the 


Figure  14.1.  Angular  variation  of  the  inner  Nusselt  number  (f  NurnKin 
for  G  =  .200  and  Pr  =  .706  (Nurmin  =  5.4836)  LUINU 


first-order  upwind  method  was  not  able  to  resolve  the  multicellular 
flow  field,  as  indicated  in  Figure  14.2.  This  was  most  likely  due  to 
artificial  viscosity  effects  associated  with  first-order  upwind- 
differencing  schemes.  Thus,  it  appears  that  at  least  second-order 
accuracy  is  needed  to  capture  multicellular  flow  transition. 

14.2.  Mesh  Resolution  Study 

This  study  examines  three  different  gap  widths  with  a  mesh  of 
31  x  102  nodes  for  both  the  complete  and  half-annulus.  Thus,  the 
number  of  angular  nodes  in  the  half-annulus  is  approximately  doubled 
compared  to  that  in  the  full -annul us. 

For  the  half-annulus,  vertical  center-line  symmetry  was  assumed. 

This  assumption  was  proven  valid  by  the  symmetric  multicellular  flow 
patterns  shown  in  Figures  6.3c  and  6.5  for  G  =  .200  and  .100, 
respectively.  In  Figure  14.3,  the  angular  variation  of  the  inner-wall 
vorticity  is  plotted  for  the  full  and  half-annulus  at  Ra  =  2,500,000 
and  G  =  .100.  The  results  of  the  half-annulus  mesh  are  nearly 
identical  to  those  of  the  full -annul us,  within  one  percent  relative 
error  for  this  pretransitional  case.  These  results  provide  confidence 
that  the  half-annular  calculations  are  indeed  valid. 

Considering  the  above  symmetry  discussion,  the  following  study 
was  performed  (note:  the  convergence  constraint  was  1  x  10"^  for  all 
test  cases  shown  below  in  Table  14.1). 

Based  on  the  results  of  this  table,  it  appears  that  the  numerically 
predicted  transitional  Rayleigh  numbers  (associated  with  the  full-annulus 


sselt  number  (7  NurnMn)  at  Ra  =  1,000,000 


Figure  14.3.  Angular  variation  of  the  inner-wall  vorticity  at  Ra  =  2,500,000  for 
G  =  .100  arid  Pr  =  .706 


Table  14.1.  Mesh  resolution  comparison 


G 

Mesh 

size 

Geometry 

RaTR 

Rah  _ 
b-a 

Percent 

relative 

error 

.100 

31 

X 

102 

full -annul  us 

2,840,000 

2840 

}  4.4 

.100 

31 

X 

102 

half-annulus 

2,720,000 

2720 

.150 

31 

X 

102 

full -annul  us 

810,000 

2730 

}  7.7 

.150 

31 

X 

102 

half-annulus 

751 ,000 

2535 

.200 

31 

X 

102 

full-annulus 

350,000 

2800 

} 1 4 . 75 

.200 

31 

X 

102 

half-annulus 

305,000 

2440 

mesh)  are  within  5-15  percent  uncertainty.  This  uncertainty  tends  to 
decrease  with  decreasing  gap  size.  In  addition,  the  transitional 
Rayleigh  numbers  predicted  with  the  half-annulus  mesh  (used  mostly 
in  generating  the  small-gap  stability  curve  of  Figure  6.12)  are 
probably  within  five  percent  uncertainty. 

It  should  also  be  noted  that  when  a  21  x  48  full-annulus  mesh 
was  used  for  G  =  .100  and  .200,  the  transition  to  multicells  did  not 
occur.  Therefore,  the  ability  to  capture  multicellular  flow  transition 
appears  to  be  very  sensitive  to  first-order  upwind-differencing  and 
mesh  coarseness. 


14.3.  Mean  Nusselt  Number  Correlations 

In  this  section,  two  correlations  for  the  mean  Nusselt  number  are 
compared  to  this  study's  pretransitional  numerical  results  for  G  =  .200 
and  Pr  =  .706,  with  a  full-annulus  mesh  of  31  x  102  node.  The 
empirical  correlations  used  are  those  of  Raithby  and  Hollands  (1975) 
and  Kuehn  and  Goldstein  (1980b),  outlined  in  Eqs.  (2.14)  and  (2.19), 
respectively. 

In  Figure  14.4,  the  above  two  correlations  were  used  to  obtain  mean 
Nusselt  numbers  within  2.5  percent  relative  error  of  those  obtained 
numerically  -  for  a  range  of  Rayleigh  numbers  between  275,000  and 
350,000.  This  comparison  is  fairly  good,  considering  these 
correlations  are  fitted  to  within  10-15  percent  uncertainty  of  the 
experimental  data. 

A  much  more  convincing  measure  of  the  accuracy  of  the  numerical 
method  was  given  in  Section  6.1.  There,  it  was  shown  that  2-D  Navier- 
Stokes  pretransitional  numerical  results  (of  local  flow  variables) 
converged  to  related  analytical  perturbation  results  in  the  limit  as 
G  -*■  0  and  Ra  -*■  °°.  Note  that  local  flow  variable  data  for  narrow  gap 
widths  could  not  be  found  in  the  literature,  hence  reinforcing  the 
importance  of  a  valid  analytical  solution  comparison. 


Figure  14.4.  Mean  Nusselt  number  variation  as 
6  =  .200  and  Pr  =  .706 


